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I. INTRODUCTION 


The growth of ocean acoustics as a theoretical and 
experimental science can be traced to a fundamental change 
in the conduct of warfare. During World War II, the 
belligerents employed scientific technology for the conduct 
of war toa degree theretofore unknown. In the United 
States scientific research was mobilized on a_e scale 
comparable to the mobilization of manpower and industry. 
Whole fields of scientific research were expanded to produce 
a technology for the conduct of the war. A portion of this 
research, the study of propagation of sound in the sea, was 
a response to a particularly effective (and technical) eneny 
weapon - the U-boat. 

The first wartime studies of underwater sound 
propagation hinted at sound transmission over long ranges 
Within a natural oceanic waveguide known as the SOFAR 
channel. A series of experiments was conducted in which 
explosive signals were propagated in the SOFAR channel over 
distances of up to 900 nautical miles (nm) [1]. In 
Subseguent experiments explosive signals were propagated 
Over paths of up to 2,300 nn j2j. Since World War II 
several uses for SOFAR propagation have been devised, 
including the location of missile impacts within wide 
oceanic areas and the study of the low-frequency ambient 
noise field. 

To predict and aid in understanding SOFAR propagation, a 
number of theoretical nodels for the transmission of sound 
in the ocean have been developed. The theoretical 
treatments have fallen generally into two broad catagories: 
ray-tracing and nornal mode techniques. Ray-tracing 


techniques are based on the Eikonal equation, 
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2 2 
Tere N= fc se(z) ] , (1.1) 
QO 


where A is the ray function, N an index of refraction, c is 
O 


a reference sound speed, and c(z) sound speed as a function 
of depth. The Eikonal equation is tTestricted to high 
frequencies in which the sound speed fluctuation is small 
over one wavelength. Normal mode technigues, more suited to 
low frequency sound, are based upon the depth-dependent 


Helmholtz equation, 


2 2 2 
CV +w /o(z) J] ¥ = 0, (1.2) 


where ¥ ais the acoustic velocity potential and w is the 


angular frequency of the acoustic disturbance. 


A. DEVELOPMENT OF THE NORMAL MODE TECHNIQUE 


Since the work to be presented is based upon the normal 
mode technique, a brief overview of the development of 
underwater sound transmission models based upon normal mode 


theory is in order. 
1. Early Theory 


During world War EE Pekeris developed two 
Mathematical approaches which became the foundation for much 
of the iater work. In 1986 Pekeris published a treatment of 
sound propagation ina layer characterized by a variety of 
sound speed gradients {3}. Although he was primarily 
concerned with sound penetration into the shadow zone caused 
by a negative sound speed gradient, he devleoped asymptotic 
solutions to the wave eguation for a number of sound speed 


Gradients. For the case of a linear sound speed gradient, 


c=c (1 + aZ), Cle 3) 
O 
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he derived a solution which reduced to a Fourier-Bessel 
Series, which converged rapidly in the shadow zone and 
slowly in the insonified zone. In addition, he investigated 
a number of variable sound velocity profiles which included, 


2 2 
c =c (1 - az) , (1.4) 
Oo 


for which the solution is a combination of Hankel functions 
of order one-thitd. 

In 1948 Pekeris published a method of determining 
the normal modes and resultant sound field in both 
two-and-three-layer isovelocity media [4]. To evaluate the 
sound field he applied a Hankel transform in the complex 
horizontal wave number plane to a vertical wave function JU, 
which resulted in an expression for the acoustic velocity 
potential 

y = fos x ak. (1.5) 
Ou. wr 
The wave function displayed a number of singularities and a 
branch cut. The solution consisted of a Summation of normal 
modes, represented by a series of simple poles, and a ground 
wave in the lowest layer, represented by a branch cut (see 
Fig 1). The branch cut extended from the branch point 


corresponding to oe the sound speed in the bottom, toi . 


This particular cut necessitated the evaluation of 
additional complex-valued poles between the branch cut and 
the positive imaginary axis (the so called "leaky modes"). 
The "leaky" modes were seen to be important only at short 
ranges. Pekeris evaluated the integral by the calculus of 
residues for the poles and the evaluation of a branch cut 
associated with the sound speed in the lowest layer. 

There were two Significant modifications to the 
Pekeris technigue which soon followed. Tolstoy reformulated 


the Characteristic equation for the modal eigenvalues 
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FIGURE 1 - Path of Bessel Integrals in Complex Plane 
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{location of the poles in the complex wave number plane) in 
terms of reflection coefficients for both upgoing and 
downgoing waves [5]. Tolstoy considered the effects of 
horizontally stratified layers above and below a depth of 


interest, Zz. A wave travelling up past 2 will be 
O O 


reflected with a reflection coefficient Rf. Similarly a 

plane wave travelling downward past 2z will have a 
O 

reflection coefficient R&R. Tolstoy showed that for an 


undamped mode representing the propagation of energy in the 


Waveguide, the reflection coefficient above and below z 
O 


must be related by 


EA RY = 1. (1.6) 


R# and Ry are functions of the sound speed profile anda 
horizontal parameter (such as the horizontal wave number or 
the horizontal phase speed). Since Eq. (1.6) is satisfied 
only for certain discrete values of this parameter, it 
represents a characteristic equation for the nodal 
eigenvalues. This general characteristic equation enabled 
normal mode analysis to be applied to a wide variety of 
environmental conditions. 

The similarity between the normal mode problem in 
acoustics and the SchrSdinger equation in quantum mechanics 
WwaS recognized and provided a tool for further development. 
Bremmer published an analysis of the Wentzel-Kramer- 
Brillioun (WKB) method of quantum mechamics aS an 
approximation to wave propagation in an inhomogeneous 
medium, opening the way for the use of the WKB solutions and 
characteristic equations as soliutions to the vertical wave 
equations [6]. 

The second change to the Pekeris technigue came when 
Ewing, sJardetzky and Press [7] suggested an alternate path 


for the branch cut (Fig 1). The Ewing, Jardetzky and Press 
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branch cut extended from the same branch point, to the 
origin, and then along the imaginary axis. This eliminated 
the need for evaluation of the complex-valued poles, as they 
were shown to be displaced to an adjacent Riemann sheet. 

The Pekeris model and its descendents were used in 
many theoretical applications, including shallow water 
propagation over fluid and elastic bottoms [8 - 15]; and was 
also used to study the effects of attenuating [16], sloping 
{17}, and rough bottoms [18] on mode propagation. 

Up to this point a major limitation for normal mode 
and many ray~tracing techniques was the required assumption 
that sound speed be a function only of depth (to separate 
the range and depth-dependent equations). In the ocean the 
sound speed profile varies gradually (and at times abruptly) 
With range. The applicability of the normal mode technique 
and, to a lesser extent, the ray~-tracing technigue over a 
horizontally varying medium was broadened by Weston in 1958 
with the development of mode and ray invariants [19]. 
Weston observed that a mode may be assumed to keep its 
identity for slow horizontal changes in the sound speed 
profile. This is known as the adiabatic assumption (after 
its analogue in time-dependent perturbation theory of 
quantum mechanics), implying that the energy propagated by a 
particular mode remains associated with that mode. A more 
recent paper by Milder expands upon the concepts of mode and 
ray invariants and derives a test for the applicability of 
the adiabatic assumption [ 20}. The adiabatic assumption 
allows local horizontal stratification in an environment 
which changes slowly in the horizontal. Thus the normal 
mode method now became applicable to a wider variety of 
Situations provided the horizontal oceanic variation was 
gradual over the long ranges. In addition, although not 
Strictly reguired, the assumption of local horizontal 
stratification has been found useful in some ray-tracing 
work [217]. 
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In the late 1950's and early 1960's there began to 
appear a number of numerical normal mode programs which took 
advantage of the digital computer. Many of the nodels were 
based at least in part on known solutions for certain 
classes of sound speed profiles. 

In 1959 Tolstoy and May developed a normal mode 
computer model based upon a series of horizontal layers in 
which the inverse of the sguare of the sound speed varied 
linearly with depth 


Z 2 
c = c (1-az) A (1.7) 
O 


In such layers Tolstoy and May approximated the node 
solutions (Bessel functions of order one-third) with 
twelfth-order polynomials [22]. For a given horizontal wave 
humber they started with both a surface and bottom boundary 
condition, then iterated from the boundaries to some 
interior depth. The Bessel functions were matched between 
layers by invoking continuity of acoustic pressure and 
vertical particle velocity. At the interior depth the 
slopes of the solutions iterated from the surface and from 
the bottom were compared. Only certain values of the 
horizontal wave number, which defined the eigenvalues, gave 
first order continuity across the interior depth. The 
solutions thus formed defined the normal modes. 

Peter Hirsch studied axial sound focusing in the 
SOFAR channel using an analytic representation for the sound 


speed profile [23]. He assumed a sound speed profile 


parabolic inc and symetric about the channel axis, 


2 Z 2 ee 
Cap eeu (i=a 2.) >>. (7.8) 
0 


For a pulse source this representation of the SOFAR channel 


Gave successive focusing at about 6 nm intervals. Williams 
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and Horne did another study of SOFAR channel axial focusing 
uSing a polynomial representation for the sound. speed 
profile [24], 


Z Z Z 3 Gu =-7 
c =c @ieaz=anz =Cz -~daZ i} ° (1.9) 
oO 


This polynomial representation of the sound speed profile 
gave a better fit to the actual SOPAR channel, allowing 
representation of its asymmetrical character. Solving for 
the eigenvalues and normal modes, they described the sound 
field for a continuous wave (CW) source. Williams and Horne 
found a much smaller degree of axial focusing than Hirsch; 
they also determined that focusing effects were sensitive to 
the detailed structure of the sound speed profile. 
Deavenport studied axial focusing uSing an Epstein profile, 

= 2 

c = A sech (z/h) + B tanh(z/h) + D, (1.10) 
with solutions which are hypergeometric functions [25]. He 
found a greater degree of focusing than Williams and Horne, 
but less than from Hirsch's symmetric profile. 

Bucker and Morris studied sound transmission in the 
surface duct, also using an Epstein profile for the sound 
speed { 26 }. The results compared favorably with 
transmission characteristics derived from experiment and 
Cay-~tracing models. 

In 1970 Williams published an outline of a normal 
mode finite difference technique [27] which started fron 
either one or both boundaries. Using a finite difference 
equation based upon the vertical Helmholtz equation, the 
method iterated across the water column to either the other 
boundary or an internal depth. An iterative approach was 
used to find the mode eigenvalues by satisfying either the 
Opposite boundary condition, or continuity of acoustic 
pressure and vertical velocity at the internal depth. A 


number of other models used this technigue, - notably the 
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model of Newman and Ingenito (28], and the model by Kanabis 
{29]. The Newman and Ingenito model started with a fluid 
boundary condition at the bottom and used a finite 
difference technigue to compute the eigenfunction across the 
water column to the surface. [The eigenvalues were searched 
by finding those values of the horizontal wave number which 
most nearly produced an eigenfunction of zero at the 
surface. The Kanabis model was similar except that it used 
a higher order finite difference equation and employed the 


adiabatic assumption for horizontal variation. 


In May 1973 the Acoustic Environmental Support 
Detatchment (AESD) of the Office of Naval Research held a 
workshop on non~-ray~tracing acoustic propagation technigues 
{ 30}. The results from seventeen acoustic models were 
compared for calculations using a variety of sound speed 
profiles. A synopsis of the models and their results is 
given by Ref. [30]. A distinction was made between purely 
normal mode {in the sense of a complete set of orthonornal 
functions) and residue series models (where the modes are 
derived from the poles of the wave function). At the AESD 
workshop a new model which represents a distinct departure 
fron previous methods was introduced. The parabolic 
equation method, outlined by Fock for electro-magnetic 
propagation in the atmosphere [31], assumes the wave 
function to be a space dependent function, V, multiplied by 
a Hankel function. The assumed soiution is applied to the 
Helmholtz equation, and terms in the second derivative of JV 
With respect to range (r) are dropped (this assumes that the 
Waves are predominantly radial - that is the corresponding 
FayS are at relatively shallow angles). The following 


parabolic differential equation results: 


a 
idV + d2yv + N (c,z)-11¥ = 0, 1.11 
se * sy * HN e201) a! 
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where k is an average wave number and N is an index of 
refraction. Since this is a parabolic differential 
equation, with a given starting solution, and surface and 
bottom boundary conditions, the solution can be solved in 
Successive range steps. The method allows full range 


dependence for the sound speed profile. 


B. OPERATIONAL MODELS 


The UG. S.~ Navy has a need to calculate low-frequency 
propagation loss over long ranges (in excess of 1000 km). 
The Navy currently has variations of four standard 


operational transmission loss models and a number of 


experimental models. The operational models include the 
following: 
a FACT - The Fast Asymptotic Coherent Transmission 


Model uses modified ray theory and is primarily for 
passive sensor applications at low to intermediate 
frequencies (50 to 2000 Hz). The FACT model computes 
propagation loss out to ranges representing several 
convergence zones [32]. FACT has been extended in an 
experimental version with multiple profiles and a 
bathymetry composed of linear segments along the line of 
bearing. The ray invariants of D. M. Milder [20] are 
used to connect rays between profiles. 

a NISSIM - The Navy [Interin Standard Surface Ship 
Model is a ray theory model which is used primarily for 
Short ranges and active sensors [33,34]. It nas had some 
limited application to lower frequency passive sensors. 

s RP-70 — This model is a multiple profile ray-tracing 
model which uses linear sound speed gradients in both the 
vertical and horizontal directions. It is used for 
passive sensors and great ranges. 


a APE - The Acoustic Parabolic Eguation model, based 


fe 





on the parabolic equation method (Eg. (1.10)), has 
recently become an operational aodel. It iS primarily 
used for very low frequencies and highly range dependent 


environments. 


C. SCOPE OF THIS WORK 


1. Objectives 





The objective of this work is the creation of a 
computer model to estimate the transmission loss of low 
frequency CW sound over long ranges and through varying 
sound speed profiles. The model is intended to be 
relatively fast, sacrificing exactness (in the mathematical 
sense) for computational speed. The underlying theory is 
based heavily upon the WKB (Wentzel-Kramer-Brillioun) 
approximation with the eigenvalues, eigenfunctions, and 
normalization all based upon the WKB formulae. The model is 


intended for the following tasks: 


a Calculation of the detailed (fully coherent) 
transmission loss for ranges up to a few hundred nautical 
Miles. In this case any sound speed variation with range 
must be slow enough t9 permit use of the adiabatic 
assumption. 

a Calculation of the incoherent Or averaged 
transmission loss for long ranges (thousands of nautical 


Miles). 
2. Departure from Previous Works 


The following are a number of the major differences 


between this work and previous normal mode models: 


e The sound speed profile is assumed to consist of 


layers {up to 50) in which the sound speed varies linearly 
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in Ee (Eq. (1.6)). Previous models have used the analytic 
solution for such layers to form the eigenfunctions 
(modes) and find the eigenvalues (Tolstoy and May [22], 
Gordon and Petersen [30], and Bartberger and Ackler [35]). 


2 
In this model the c linear profile 1S primarily used for 


the ease with which the vertical wave number may be 
integrated - an operation required for many of the WKB 
calculations. The analytic solution (the Airy functions) 
is used as a mode solution only where the WKB 
approximation breaks down. 

»s The sound speed profile is divided into separate 
ducts corresponding to different families of modes. The 
calculation of eigenvalues and modes is carried out family 
by family. 

e Within each family the behavior of the eigenvalue as 
a function of mode number is approximated by a fifth-order 
polynomial. This polynomial provides the first estimate 
of the eigenvalues for each mode. For calculation of the 
averaged or incoherent transmission loss at great frange, 
this estimate is sufficiently accurate for computation of 
the eigenfunctions. For detailed transmission loss at 
closer range, the eigenvalue is refined by a maximum of 
five iterations. 

s The WKB approximation has been extended to 
incorporate the effects of a.barrier between multiple 
sound channels, and the effects of the surface and botton 
boundaries. The eigenfunction is found by a combination 
Of WKB, Airy and extrapolated functions. 

s A range dependent environment has been included by 
the use of the adiabatic assumption. A previous 
difficulty with the adiabatic assumption for normal node 
use has been the ambiguities for modes in different sound 
Channels of multiple channel profiles. TLasoweG tf 2oule y 


has been overcome with the use of a formalism which 
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matches modes within equivalent channels and approximates 


the transition of modes between channels. 
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II. THEORY: HOMOGENEOUS WAVE EQUATION 


—_— aa ae _— ep ee ee ee ee ee aS 


In what follows it is assumed that the ocean ls a 
horizontally stratified medium of constant density, with a 
perfectly flat and pressure release surface, and a flat 
bottom. The input data for most ocean acoustic models 
represent a linearly-segmented sound speed profile. fhe 
linear segments do not represent the detailed structure of 
the ocean sound speed profile; they are a convenience for 
interpolating between the data. In ray trace models this 
linear gradient takes on added significance because of the 
Simplified ray path thus formed. Because of the analytic 
solutions available, when using normal mode techniques it is 


convenient to approximate the sound speed structure between 


data points as linear in c . The profile will be assumed 


to consist of layers within which the sound speed varies 


-1/2 
ae ce tt = a (z-z)j) ,2z. <z2<2z,, (2a 
1 1 i 1+1 


where c andz are the sound speed and depth at the top of 
1 a: 


th 
the 1 layer (Fig 2). If w is the angular frequency of the 


acoustic disturbance, then the sound speed profile can be 


represented in terms of the wave number, k. For the profile 
2 
of Eq. (2.1) k varies linearly with depth, and the gradient 


2 
of k , y, is defined by 


2 Zee 2 
Fee oeeyGeea=ek = y {Z7>Z ), Z < 2S Z. Cf (2.2) 
Be 1 a 1+1 
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FIGURE 2 - Sound Speed Profile 
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The horizontal wave number k and the vertical wave number 
k are related by 
Z 


k =k <-k. (2a) 


The vertical wave number is real for k > k , andis 


be 
imaginary for for k<k ; 
E 
k = -i 6, 
Z 
where 
Ze 2 2 
B =k -k (2.4) 
E 
moc K < ki. 
c 


The theoretical development which follows can be divided 
into five steps outlined below: 


a The homogéneous wave equation is established, and 
from it 1s separated a depth-dependent equation. The 
depth-~dependent eguation must be solved in terms of 
eigenvalues and corresponding eigenfunctions {the normal 
modes). 

a TwO approximate solutions’ to the homogeneous 
depth-dependent equation, the WKB solutions, are derived. 
One oscillatory solution corresponds to the trapped field 
in the waveguide, the other, exponentially decaying, 
represents the diffraction of sound at the edge of the 
waveguide. The WKB solutions fail where they approach 
each other at the edge of the channel (turning points). 
From one of the WKB solutions a characteristic equation is 
developed, whose solutions are the modal eigenvalues. 

a A third approximate solution to the homogeneous 


eguation is derived. This solution is used to connect the 
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two WKB solutions. AS a consequence, for each eigenvalue 
an unnormalized eigenfunction can be determined as a 
function of depth. 

«a The treatment of the homogeneous wave eguation is 
completed with consideration of the effects of the surface 
and bottom boundaries. 

# The final step is the introduction of an 
inhomogeneous term (representing a point source) into the 
wave equation. The inhomogeneous equation 1s then solved 
by the residue series of a Green's function which is 
formed from the previously developed homogeneous 
solutions. Solution of the inhomogeneous eguation allows 
determination of the mode amplitude to match a point 


source. 


A. HOMOGENEOUS WAVE EQUATION 
The homogeneous lossless wave eguation for velocity 
potential Y 1s 
y2y 2-1 d2yY¥ = Q, (2.5) 
The acoustic pressure, P, is related to the velocity 
potential by 


P= — 7e1dv a. 
dt 


For a monofrequency acoustical disturbance with time 


dependence exp (-1i wt), 
P= iwo¥. 
The boundary condition at the surface requires that 
z= 0) —0r (2.6) 


where zis the depth below the surface. The boundary 
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condition at the botton ee Gane wWertten ih the £fornr 


d =u, i 
it a { ) 


rE (a 


where , will be determined by the particular bottom model 
used. In addition, certain restrictions must be placed on 
yi s 
ee | must represent an outgoing wave in the 
horizontal direction and an exponentially decaying wave at 
great depth, so that { ¥ | *+ 0 as r*>?% orzr>e.,. 
2. There must be no discontinuity in pressure and 
particle velocity in the water column so that pY¥ and the 


gradient of ¥ must be continuous. 


i. fhe Helmholtz Equation 
If the time dependence of y is explicitly stated, 
¥ = 6 exp (-iwt), | (2.8) 


where 6 ais a function only of space, then substitution of 
meme (2.8) into Eq. (2.5) yields 


V2 6 + »o = Q. 2a 

a (2.9) 
Expanded in cylindrical coordinates with the origin at the 
sea surface and z positive downward, Eq. (2.9) becomes 


teom(Gdd) Sede t + ure = 0, (2.10) 
padre ear dz< ce 


for cylindrically symmetric motion. 
2. Zhe Hankel fransfora 


EES SSS SS a ee 


One wavy to separate the depth dependent terms in 
Eq. (2.10) is to apply a Hankel transforn. The Hankel 


transform variable, k , represents a horizontal wave number, 
c 
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and transforms the solution from a depth and range space 
(z,c) to a solution in depth and horizontal wave number 


Space (z,k ). The transform pair is defined 
c 
U(z,k ) = fe (2b) vo (kor) EF dar, (25,11) 
E oe « 


and 


d (z,L) - | U(z,k ) 3 (k rc) k dk. (2.12) 
Cc Oo Cc Cr bo 


Application of this transform to Eg. (2.10) yields 


Z 
d2U + 2-k O(a ©) Or 2213) 
dz2 Lue rc NCU 1s 
This is a homogeneous differential equation with the 
following conditions imposed as a consequence of the 
conditions on y : 


a) U(0O,k ) 0, the free surface (2.14) 
a 


boundary condition. 


b) 1 dU = un, at the bottom, where z=2z. 
U dz b 


c) AS z+o , U must represent a decaying wave. 
qd) o U is continuous. 


e) dU is continuous. 
dz 


ee Natur 


of the Solution 
The function U(z,k ) defines one normal node, which 
iS 
Cepresents the depth component of an acoustic wave 


propagating outward With cylindrical divergence. The 


horizontal wave number defines the phase speed, 


res =yw / k a (2515) 
Dp Er 


The phase speed represents the speed with which surfaces of 


constant phase are propagated in range. 
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Inspection of Eq. (2.13) reveals some gualitative 
properties of the solution. For k> k the vertical wave 
r 
number, k 1s real and the solution is expected to be 
Z 


oscillatory in nature. Por k < a the vertical wave number 
is imaginary, and the solution is expected to be exponential 
in nature. Regions for which 5 is real (oe and which 
are bounded by depths with imaginary ie are) form an 
acoustic channel or waveguide. Within such a channel the 
acoustic waves are refracted between depths at which s = 0, 
or reflected from boundaries. 


The ocean normally forms an acoustic waveguide for 


phase speeds less than the speed of sound in the bottoan, . 


At these phase speeds the modes correspond to waves which 
are either completely refracted or strike the bottom at a 
graZing angle less than the critical angle. For the 


corresponding values of k (k ea Eq. (2.13) with the 
ig rc 


conditions of (2.14) forms a Sturm-Liouville boundary value 

problem over the interval z = 0 to z= os AS a consequence 

there is a finite number of eigenvalues for the 
’ 

parameter c and corresponding eigenfunctions, ve which 

solve the Sturm-Liouville problem [1]. Thus the range of k 


> S is termed the discrete spectrun. 


At phase speeds greater than the sound speed in the 
bottom, the differential equation and boundary conditions 


are met for all phase speeds or k . For these high phase 
i 


speeds there 1s a continuous and infinite set of solutions 
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Mor ail k <k. 
Cc b 


The reader familiar with the Schr&Sedinger equation 
of quantum mechanics wiil recognize that the sound velocity 
profile is an analogue to a potential well. The sound speed 
profile itseif corresponds to the potential function, the 
phase speeds correspond to the energy levels, the surface 
boundary to an infinitely high “wall", and the bottom 
boundary to a "wall with height that corresponds to the 
sound speed in the bottom [2] [3]. 


B. THE WKB METHOD 


In this section the WKB method, an approximate solution 
originally developed in quantum mechanics, is applied to 


solution of the wave equation for normal modes in the ocean. 
1. The WKB Solution 


RE —p-mem =a 


To derive the WKB formulae a method used by Schiff 
is adapted to oceanic parameters [4]. By defining a scale 


parameter, 
h= if, 


and an index of refraction, 
2 2 z 
Nese = al7c »l, 
p 
the vertical wave number can be re-written, 
k = N /h. (210) 


Assume that the depth function has the forn 


ie exp{ iS {z) J. (2.17) 


35 





where S$ is some (Unknown) function of depth. Substitution 
Seeeeas. (2.16) and (2.17) into Eq. (2.13), gives 


Z Z 
ihs'' - st +N =0. (2.18) 


Expand S in powers of h, 
Se-woeet ASS ff fae 3 (2.19) 


and equate terms in the first two powers of h, resulting in 


two equations 
a +N = 0, (2.20) 


and 


eS ee 2 So! a) 
0 Ca 


Integration of the first equation yields 


Zz 
S {z) = £ f Naz 3, (Zea) 
0 Z 9 


and substitution into the second gives 


S (z) = i log N. DED 
ae Same? ( ) 


Substitution into Eq. (2.17) gives the two WKB solutions. 
Por k > k , 
r 


/2 Z Z 
U(z) = k {a cos (| k dz) + b sin(| Zita (2623) 
yA oO 2h Zz oO 2A Z 


This solution will be referred to as the WKB oscillatory 
solution and can also be written 
-1/2 ze 
U (z) = k eos (7) kK-dz-0 ), (2.24) 
I Z ae: O 
0 


where 
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@ = arctan (b /a). (2525) 
O Oo Oo 


For k < k the WKB exponential solution is obtained, 
i 


Z 
ee 
= {c exp( /8dz) +d exp(-[ 8 42)}. (2.26) 
° Z 6 . Zo 


2. WKB Validity Requirements 


— =< éo§ Saag =e EE a = = 


In order to expand S in powers of h (Eq. (2.19)) 
successfully and to separate the first two terms, the second 


tern, aS must be small in comparison to the first, S. 


This requirement is met if [5] 


[hse /s*| = Idi d_k Gar. (2a 


A WKB validity factor, FPF , will be defined as equal to the 
W 
two-thirds power of the expression in Eq. (2.27). (The 


two-thirds power is used to facilitate comparison with the 
Airy validity factor derived in a later section.) From 
Ege (2.27) with the help of Eq. (2.2) 


273 ine 17-3 


F = zig | Y Gri; (2223) 


1 

W q 

where z is the turning point depth. Thus |{F | must be much 
ite W 

less than one for the WKB solution to be ae_=e good 

representation. In terms of the horizontal wave number the 


expression for F becomes, 
Ww 


2/3 22 ey 
PF = 14 Veen | hn ene (2.29) 
W q Lc 


ay 





of the WKB Solutions 


3. Properti 


if 


The solution to Eq. (2.13) is characterized by an 
oscillatory <tunction (Eg. (2.24)) at those depths for which 
k 1s real. This function corresponds to the summation of 

Zz , 
an upgoing and a downgoing vertical wave. At those depths 
for which k 1S imaginary, the solution represents the 
Z 
lateral wave in a shadow zone. For the reference depth z 
oO 


in Ba. (2.24) and (2.26), the turning points or a 


reflective boundary will be used. 


The @ of Eq. (2.24) cepresents a phase change as 
O 
the vertical wave is refracted at the turning point or 


reflected at a boundary. Considering the oscillatory 
solution as the combination of an upgoing and a downgoing 


wave, 98 represents the phase change for only one of the 
O : 
waves. Thus the total phase change experienced by a wave as 


it travels through a turning point or is reflected from a 


boundary is twice 9. 
oO 
At the turning point F becomes singular; thus there 
W 


1s some region bounding the turning points for which the WKB 


solutions are not valid. 


4. WKB Characteristic Equation 





The homogeneous differential equation (Eq. (2.13)) 
can be solved and both boundary conditions satisfied only 
for certain eigenvalues of k . After developing the WKB 

B 


approximations, the next step is to use the WKB oscillatory 
solution to determine a characteristic equation for the mode 


eigenvalues, k 3 
Berea 


38 





Consider some arbitrary depth, z , between the two 
1 


turning points for a specific k (Fig 3). Assume 
En 
U(z,k ) is a solution to Eq. (2.13) which neets all the 
eit 
conditions of (2.14). From Eq. (2.24), at depths above oe 
the solution must be 
Z 
-1/2 : 
U (z ~) = k cos ( [x adz-9 ), 
ian 4 Z Z u 

tu 

where Z is the upper turning point depth and 9 is the WKB 
tu u 


phase angie at the upper turning point. Conversely, if the 


lower turning point ot 1s used as the reference level, the 


solution below ee must be 


mel 
-1/2 
Ue(2 +) = k cos(/]k dz-9 ). 
2 1 Zz Zz 1 
al 
For the two expressions to be eguivalent, they must be 
continuous through the first derivative at Zie This 
requirement is met only if 
Zz 
tel 
[x dz-9 -¥ = nw, i= Oe vecigters s (2.30) 
; Zz u Ll 
tu 


This is the characteristic equation for the WKB 


approximation. If 98 and ot equal7w/4& this characteristic 
u 


eguation becomes the Bohr-Sommerfeld quantization rule for a 


potential well [6]. 


Consider a wave which prepagates upward from z_ to 


the upper turning foint ae , is refracted downward to ro 
u 


and then is refracted upward back to depth a The 


a 
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FIGURE 3 - The WKB Solutions 
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characteristic equation implies this wave has undergone 


phase change which is an exact multiple of 2 . 
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III. THEORY: AIRY SOLUTION AND CONNECTION FORMULAE 


ee oe ee 





After deriving the WKB solutions, which are good 
approximations away from the turning points, the next step 
is to derive a solution valid at the turning points. This 
solution, called the Airy solution, will be used to match 
coefficients between the two sets of WKB eigenfunctions 
determined from Eg. (2.23) and (2.26). The derivation will 
follow the method used by Lauvstad [1]. 


A. THE AIRY SOLUTION 


2 
For some region about the turning point, k varies 
r 


linearly with depth and Eg. (2.13) may be expressed 


d2yU - y(z-z )U = 0 (3.1) 
azz t ‘ 
where y is defined in Eq. (2.2) and ze is the turning 


point. With the introduction of a new variable, 


1/3 
x = y Zee (3-2) 


the second derivative with respect to depth becomes 


2773 
Zee a (dee x) =) ¥ qe, 
Ze dz dxdz ax= 
resulting in the Airy equation, 
d2e—sxi = 0 (3.3) 





Solutions to this differential equation are Bessel functions 
of order one-third. The solution to Eg. (3.3) may be 
expressed in terms of a particular combination of these 
Bessel functions, the Airy functions Ai and Bi [2] 


U(x) = A Ai(x) + B Bi(x). (3.4) 

O O 
With approximate expressions for the Airy functions for 
large and small argument, the two WKB solutions may be 


matched to the Airy functions and, in the process, matched 
across the turning point. 


1. Airy Solution for Small Argument 


For small x the expressions for Ai and Bi can be 
represented as ascending series [3]. Taking the first two 


terms of those series, we have for |x] << 1, 


U ene ana B c B-A) X, 3.5 
pen ee x + piv ‘ >: X ( ) 


where 
cr = 0.35502805, s = 0.25881940 . 
This expression, the Airy solution, 1S a good approximation 


at the turning point, precisely where the WKB solutions 
break down. 


Meee Alry Validity Factor 


The WKB solutions are good approximations at some 
depth beyond the turning points. In addition, there is sone 


interval about the turning point cm Within which the Airy 
Solution 1S a good approximation. The requirement that {x{] 


<< Tin Eq. (3.5) may be imposed through the use of an Airy 
maenaity factor, F : 
a 
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ry 
it 


ese. 
|x| 1 oy az) << 1, (3.6) 


or 


-2/3 2 Z 
{x y | k - k fe (3-7) 
a Cc 


F 


The validity factors F and F are plotted as functions of a 
a W 
1/3 Has 
scaled depth, Y Meme (2g 7h). Y£ limits of 0.3 and 


0.448 (0.3 to the two-thirds power) are used for F and F 
a W 


respectively, the Airy solution is a good approximation for 


x < 0.3 and the WKB solution is a good approximation for 
moe > 0.87. There is a “gap" for 0.3 < {xj < 0.87 within 
Which neither the Airy nor WKB solutions is a good 


approximation. - 


3. Airy Solution for Large Arqument 





To determine the correspondence between the Airy 
solutions and a particular WKE solution, the coefficients of 
the Airy functions for large argument and the WKB solutions 
are matched. Introduce a variable (solely for notational 
convenience), 


a2 


¢ = X : 


Z 
3 
for {|x| >> 1, the Airy solutions may be written in terms of 


the asymptotic expansions for Ai({x) and Bi(x) [4]. For the 


Shadow zone region (k> k), the Airy solution, for 
rc 


asymptotically large argument, is 


isis 41s fi'ae +Be j. (3.8) 


G4 





é 


‘ yqdeq pateos jo wot youn e se PTz9qytTAD AVTPTTVA - hf ARNOT 
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For the oscillatory region ( k < k), the Airy solution for 
c 


large argument 1s 


1/4 


U(-x) = [B cos(c+n) + A sin (sn) J. (3.9) 


ax 
& 


B. CONNECTION FORMULAE 


The next step is te equate both WKB solutions, oe and 
U , to the Airy solution, JU , then match the WKB 
iI ey 


solutions across the turning point. In addition, the case 


of a barrier (sound speed relative maximum between two 
channels) is considered, and the WKB solutions for the 
channel on one side of the barrier are matched with the 
solutions on the other side of the barrier. 


The WKB solutions are a function of the integral of the 


vertical wave number. For a single linear segment (in k ) 
this integral is expressed 


Z 
0 1f2 3/2 
ik {dz = 2 + Jz~-z | ae (3110) 
2 Z 3 oO 


For notational convenience a constant, g, is defined by 


1/6 -1/2 
( T 


gS Sy ) - 


1. WKB Oscillatory and Airy Connection 





By matching coefficients, the oscillatory WKB 


solution Se and the Airy solution a) are related by 





at A +B), B=) 1 A= Dae 3.11 
os ar ( ‘ 3! ; oy ( 3 3! ( ) 


A = _1 a+b), Be eet a-b). 
725 ( a ot ( . 3) 


0 6729 


2. WKB Exponential and Airy Connection 


The exponential WKB solution oe and the Airy 


solution (U ) are related by 
ae: 


d = g Ao, c =g 8B, (3.12) 
oO oO 


22) a ee = GEES CEE ree SS EE ee eee =] Se Gea Ga ee ee Set 


Binoy; bpyecompanang EG. (3.114) and (3.12), the 
two WKB solutions are matched by 


= 4 2a + ' b = 4 2dad- ’ 3213 
%e SL : Oo ce oO Az Oo oA : 
= 4 -b), d = 1 Dye 
a Je oe a Oo Zo oo By 


The 9 of Eq. (2.20) can now be re-defined as 
O 


@ = arctan (b /a ) = arctan [ (2d -c )/(2d tc ) }. (3.14) 
oO o 60 Oo 8600 Oreo 


C. BARRIER CONNECTION FORMULAE 
When a sound speed profile has two relative minima (two 


Channels), it is necessary to describe the lower turning 


point phase angle orf the upper channel, Sa in terms of the 
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upper turning point phase angle of the pee channel, 90, 
and S (Fig 5). Consider a wave which approaches the lower 
turning point of the upper channel, ane The wave is 
refracted and vertexes in the vicinity of the turning point, 
mrs Some of the wave energy “leaks” through the barrier and 
enters through the upper turning point, Ce into the lower 


channel. 


If the barrier 1s thick enough for the wave to be 
considered as progressing from an Airy representation to a 
WKB representation, and then back to an Airy representation, 
the connection formulae (Eg. (3.13)) may be invoked twice. 


Thus the WKB coefficients in the upper ensonified channel 
(a ,b ) in terms of the coefficients for the lower channel 
Oo Oo 


i”) Can be expressed as 


L 
a = (a -=b)e + 1(a +b )e 
( > a x! 5 . ’ 


Le =I, 
b = (a ~b )e - 1(a +b )e 3.15 
( 3 ) x! > >) ’ ( ) 


#here 2 


2 
pes GZic 


Z 
{[ 8 is the imaginary vertical wave number in the barrier, 


see Eq. (2.4) J. The WKB phase angle at the lower turning 
point of the upper ensonified channel, oe can be 


expressed in terms of the phase angle at the upper turning 


point in the lower channel, ee 
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FIGURE 5 - Double Sound Channels and Barrier 
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8 = arctan (b /a ) 
1 oe) 


L L 
= arctan tien’ ec = eco ye ey 


L L 
eae.) © + nouuscane) © }- (3.16) 


If the upper channel 1S considered the ensonified 
channel and the lower channel that into which energy is 
"leaked", the relative amplitude of the lower solution 
represents the transmission coefficient for the barrier. 
This requires that the acoustic amplitude at the top of the 
barrier exceed the amplitude at the bottom of the barrier. 
Thus, the amplitude of the WKB exponential solution at the 
top or the barrier should be larger in magnitude than the 
WKB exponential solution at the bottom of the barrier. In 
terms of the coefficients of Eq. (2.26) this may expressed 


Z iy aii 2 
fcmmtede | > ice + de jf, 
oO oO oO oO 


or 


=1, 
Fowl << ldelie 2 (3.17) 
O O 
When this is substituted into the connection formulae 
(Eq. (3.13)), a restriction is placed on the lower phase 


angle of the upper channel 


mn ~- arctanjf[c /f2d | < 8 < 7 + arctanj[c /2d |. (3.18) 
q oO oO 1 q oO oO 


Thus for a wide barrier, the WKB phase angle on the 
ensonified side of the barrier approaches "/4. As the 
barrier narrows, the phase angle may be perturbed an 
increasing amount on either side of the 7/4 value. As the 
barrier vanishes, L in Eq. (3.17) approaches zero and the 
phase angle may vary from between 18.669 to 71.34°. 
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IV. BOUNDARY CONDITIONS 


In this chapter expressions for the surface and bottom 


boundary conditions are developed in terms of the 694 OF 
O 


Ege (2.24). Also, in order to evaluate the eigenfunction 


normalization factor (derived later), expressions for the 


derivative of @ With respect to k are obtained. The 
oO E 


derivations wili be divided into three cases, based upon 


whether the sound field at the boundary is best described by 
the WKB exponential solution (Eq. (2.26)), the Airy Solution 
(Eq. (3.5)), or the WKB oscillatory solution (Eq. (2.24)). 


A. THE SURFACE BOUNDARY 


The pressure release boundary condition at the surface 
as defined by Eq. (2.14a) 1s applied to one of the three 
solution cases. The coefficients thus derived are matched 


to produce a value for @ in Eq. (2.30). 
u 


ime Case 1 


Consider a surface boundary which is in the shadow 
zone region. If Eq. (2.27) is satisfied at the surface, the 
WKB shadow zone solution (Eq. (2.26)) is a good 
representation at the surface (Fig 6). Apply the surface 
boundary condition to obtain 

-1/2 Feu mtu 
U (0) = 8 {c exp( |] gdz) + d exp(-]| 8dz)} = 0, (4.1) 
LL O ; O ; 
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rUGunreo = surface Boundary Condition - Case I 


a3 


where Br is the value of g at the surface. Using the 


notatioa 


Z 
tu 
Lu = SeaZ, (4.2) 
0 
and with the connection formulae (3.13) and (3.14), the 


surface boundary condition may be written, 


u u 
@® = arctan (b fa ) = arctan [ (2te 7 t2-e ) Je (4.3) 
u Oo Oo 

As the depth of the turning point increases, Lu increases in 


magnitude, and approaches ,/4. This result corresponds 
u 
to the 1/4 phase shift for a refracted wave in an unbounded 


medium, as derived by many authors including Brekhovskikh 
mam, schiff (2j, and Tolstoy and Clay ({3]}. Differentiating 
with respect to k obtain 

c 


a (8 eee te: a ea) 4.4 
= ~-¢se +e Us « - 
me ni ans a 


e if 


NO 
® 
Is 


Etaeicie <s | in Eg. (3.5), and the surface and 
2 
turning point occur in the same linear segment of k, then 
the Airy solution (Eq. (3.5)) can be matched to the surface 
boundary condition (Fig 7) 
c (A +/3B ) +c (/3B -A)x = 0. (4.5) 
io oO Z Oo 8060 

Using the turning point connection formulae, the upper 
turning point phase shift may be expressed, 


@ = arctan (bd /a ) 
u (oe) 
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ao 





where 
b fa ={[c (f3-1)x + c (f3+1)] / 
oOo oO Ze 1 
(c, W3+1) x tc (/3-1) ]. (4.6) 


For values of [x] < 0.3 the upper turning point phase . 
shift increases from 65.319 to 87.009 as the turning point 
approaches the surface (and the mode becomes’ surface 
reflected). This differs from the WKB exponential upper 
turning point phase shift, which varies from 45° to 72.579 
as Lu decreases from a large number to Zero. 


The derivative of 8 with respect to k is 
u is 


fy OA 
a (8) = (a tb) 2k ca (f3-1)-b (f3+1) ]. (4.7) 
qk u oO Ee 3201-0 O 
When x equals zero, which corresponds to a wave which is 


refracted with grazing incidence at the surface, 


@ = arctan [ (/3+1)/(/3-1) J, 
u 


= 57/12 radians or 759. 


This condition is analogous to the ray that grazes_ the 
surface at zero degrees, or is just tangent to the surface. 
This result can be duplicated by starting with the 
connection formulae given by Schiff [4%], rather than with 
the connection formulae (3.13). 


Since Eq. (4.5) 1s an approximation for Weak on both 


Sides of the turning point (for small values of x), 


Ege (4.6) should provide a good representation for small 
negative as well as small positive x. This corresponds to 


Waves reflected from the surface at extremely smali grazing 
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angles. 8 as a function of x increases as the denominator 
u 


of Eq. (4.6) approaches zero. This occurs for a value of 
x = Co (1-73) J / [o, (14/3) J = -0.3675516 . 


At this value, 


@ = r/2 radians = 909, 
u 
and we have the most common form of the pressure release 
boundary condition - a phase shift of 7 radians or 180 
degrees (here split equally between the upgoing and 


downgoing waves). 


3. Case III 


ie 


In this case consider a wave that is reflected from 
the surface at an angle large enough so that the WKB 
representation is. valid. Thus the boundary condition 


becomes (Fig 8) 


a 
-1/2 
U = k Cos t(y KedZ—-9 )e— 0. (4.8) 
a Z 0 Z u 
At the surface, Z = Zz and 
oO 
g9@ = T/2 or 90°. (4.9) 
u 


This is the expected result for a pressure release surface. 





The results for the three cases described above are 
. . 2 . 
Shown in Fig 9 as a function of x, for k a linear function 


of depth. The WKB validity criterion F is also displayed 
W 


(the validity criterion for the Airy solution is [x]). The 
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VSe 


Airy analysis shows that there are three distinct cases for 
the effect of the surface pressure release boundary upon the 
upper turning point phase shift. The extent of each case 


depends upon the gradient of k at the surface. Note that 
Is 
there is a “gap" between those values of |x| which are 


sufficiently small for {PF |] << 1, and those values which are 
a 
sufficiently large for |{F { << 1. This gap could be 
W 
shortened by the inclusion of more terms in the ascending 


series for the Airy functions. 
In addition, the sound speed profiles normally do 
not allow the representation of the sound speed curve 


between the surface and the upper turning point as a_ single 


2 
linear gradient in k . Usually there is a combination of 


such linear segments, which may affect the results. 


Consider a profile near the surface as a combination 
of linear segments as in Fig 10. Since the boundary 
condition at the surface must be satisfied, it must be 
determined which, if any, form of the solution may apply at 
the surface. In most cases there should be a range of 
relatively low phase speeds for which the WKB shadow zone 
solution is a good approximation at the surface. This range 
corresponds to the deep channel refracted waves (Fig 10). 
At the other end of the phase speed spectrum, there will be 
some range for which the WKB oscillatory solutioa (with a 
phase shift of 1/2) is a good approximation. In addition 
there will be some narrow band of phase speed, centered 
about the surface sound speed, for which the Airy solution 
1s a good approximation. 

The gaps between these ranges of good approximation 


could be pridged by using the full Airy function 
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FIGURE 10 - Typical Sound Speed Profile near the Surface 
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~ 


representation and then matching the Airy functions fron 
segment to segment. This is the basis of models developed 
by Tolstoy and May [5], and Bartberger and Ackler [6]. 
However, the calculation of the Airy functions for a large 
number of layers becomes a time-consuming task, and is 
beyond the aim of a fast estimate set for this model. 

If the surface boundary condition is plotted as a 


function of c , the results are Similiar to Fig 11. For low 


C the value of 6 approaches 7/4 radians or 459°. Asc is 
P ut Pp 


increased, a region is approached where the shadow zone WKB 


solution is no longer valid. As c continues to increase, 
P 


the other side of this "gap" is reached where the Airy 


series solution of case II is a good approximation. This 


case has a value of 75° for @ at the c corresponding to 
u Pp 


the surface sound speed. Case II continues to case III 

where the value of 8 is a constant 909°. In this computer 
u 

model (QMODE) the gap between the case I solution and the 


case II solution is bridged with a third order polynomial in 


r - 
E 


B. BOTTOM BOUNDARY CONDITION 


fo find the bottom boundary condition the sane 
techniques are used as for the surface boundary conditions. 
The boundary condition at the ocean bottom 1S expressed by 
Hae (2. 1%b). Assume that a value for , and the derivative 


of uw with respect to k is provided by the bottom model 
E 


used. Normally y will be a function of k . 
rc 
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As with the surface poundary condition the lower turning 
point phase change, oar Will be investigated for the three 


forms of the solution. 


1. Case I 


For modes refracted at depths well above the botton, 
the WKB exponential solution is a good approximation 
(Fig 12). The derivative (with respect to depth) of the WKB 
exponential solution (Eq. (2.26)) is 


mun Aa. LL =Li 
aq {0 ) = 8 (c e -de ) 


Oz it oO Oo 

WameCunens 446 ya (4.10) 
= ce +de 

3° O O ae 2 

where 
25 
Ll = J eaz. 
“+ 


By use of the WKB assumption in Eq. (2.27), drop the second 
term in Eq. (4.10), and obtain 


TE -LL ek -Li 
u = 8 (c e -de ) f/ {(c e +#+de ye (4.11) 
fF Oo fe) Oo Oo 


where 8 is evaluated at the bottom of the water column. 
Solve for c andd, and substitute into Eqs. (4.7) and 
oO oO 


(2.26) : 


t 9 b 2 s aa / 
an = A = -,; )e -({g +y )e 
7 Ea : [ (8 uw) (6. yu) J 


Ee 
eee re O08 nye ]- (4.12) 


Note tnat as Ll becomes larger, oh approaches 7/4 radians or 


459, which agrees with the phase shift for a ovourely 


refracted wave in a boundless’) mediun. AS Ll approaches 
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FIGURE 12 - Bottom Boundary Condition - Case I 
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zero, the turning point approaches the bottom and @ 
approaches 71.57°. 


Differentiate oF with respect to k to obtain 


i 
a9 cea Wo feet 4.13 
dk ,1 ie Ae io Oo Oo oft ( 


where primes denote partial differentiation with respect to 


k . For this case a' and b' are evaluated by 
c o ° 


LL -tLl pe 
b L1' + (2e te )B - (2e -e jut, 
oO 


a? 
O 


ee — 1 li LL -L1 
a Llt + (2e -e Bf - (2e te Piao (4.14) 
O 


oO 
it 


L1' and 8 are evaluated by 


b -1 
Ll' = @ (Ban =k g dz, (4.15) 
dk E 
742 
and 
ee) eee (4.16) 
Be ae oe ey 


The bottom model used must supply a value for u'. 
Ze Case il 
For k sufficiently close to “a the wave number at 
r 


the bottom of the water column, the Airy solution 


(Eq. (3-5)) 1s a good approximation Craeeme ls ie Fron 
Pas (3.5) 
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2/3 
Hee te (798-8) }] /{c (eae yaya + 
2 Oo Oo 1 Oo Oo 
2 
c (f3B -A )B_ ]j. (4.17) 
2 Onno: > — 


Solving for A and B_ gives 
oO oO 


2 2/3 
b fa = [(c_(Y3-1) (yg J-c Y3t1)u ¥ 1/7 
o 600 2 f 1 
2 2/3 
[o, /3*1) (YB )-c (3-1) u y je (4.18) 


For a rigid bottom y= 0, at k equal to ae a equals 1/12 
rc 
radians or 159°. fif approaches infinity (which corresponds 


to a pressure release surface), Eq. (4.18) becomes 


Eqs (4.6). For case II the derivative of e (Eq. (4. 13)) 


is evaluated by 


gm 
ay 
ii 


2 273 
-2c k (M441) y - [c_8. Wati)+c y  (f/3-1) jut, 
2. <r 2a 1 


2 2/3 
b' = -2c k (f3-1), - (¢_ 8 (/3-1)+c rae (/3+1) Jul. (4.19) 
O 2% 2 1 


3. Case Jil 





The Simplest case is when the wave is reflected from 
the bottom at grazing angles great enough to allow use of 
the WKB osciliatory solution (Fig 14). The derivative of 
the WKB oscillatory solution (Eg. (2.24)) is 


Z 
Ve o 
= k Sin( Jk dz~9 ), (4.20) 
Z Z ip 


Z 


d_ (vu 
gs‘ = 


where the terms for the derivative of k have been dropped 
Z 
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@eoein Eq. (4.10). Thus 
“ 
u = k tan {x GZ ou): 
Zz Z aE 
Z 
evaluated at the bottom, the expression for o. becomes 





@ = arctan (-u/k é 4.21 
iL op P bottom ( 
The derivative with respect to k is 
i 
a (0) Femi Fokus k /k ia? 
= + im = ° ° 
aK 1 ou Zz Z s 3 if a 





As with the surface boundary condition there are 
three ranges of phase speed for which each of the forms of 
the solution is a valid representation (Fig 15). At 
relatively low phase speeds, the shadow zone WKB solution is 
a good approximatian (case I). At phase speeds near the 


sound speed at the bottom, the Airy solution is a good 


representation (case II). At high phase _ speeds near 
cut-off, the oscillatory WKB solution 1s a good 
representation (case III). In the case of the botton 
boundary, howevel, there are two gaps between the 


representative solutions: one between case I and case II, 
and one between case II and case III. Ina method similar 
to the surface case, the gaps between these caseS can be 


bridged with two third-order polynomials ink. 
E 


C. BOUNDARY CONDITION TES? 


The eguations for the bottom boundary condition were 


tested for three test cases Dy comparing the WKB results (as 
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calculated by the program QMODE) with a finite difference 
technigue. The finite difference technique was adapted fron 
a normal mode program by the author [7]. The routine 
started with a fluid bottom boundary condition and iterated 
a solution for the depth function upward over a depth grid 
of 501 points using the finite difference equations of the 
NORMO1 program of Ref. [7]. At the grid points the sound 


_speed profile was interpolated to be linear inc . This 


was done to insure the closest possible match between 
profiles. The iteration was stopped when the lower turning 
point had been crossed and 


2 
tk (it1)-k (i) t hk 0 Oe (4.23) 
Zz 2 zZ 


where h is the grid step size. This condition allows the 
WKB form of the solution to be a fairly good approximation. 
The lower turning point phase angle was then calculated 


uSing the expression 
Z 


@ = arctan [-Ut/(k JU) j -[ dz 
1 Z z 
+1 
where U' is the derivative of U with respect to depth (this 
quantity is a part of the NORMO1 iteration scheme). The 


integral of k was calculated using the trapezoidal rule. 
Zz 
If Eq. (4.23) was not met, no comparison was made. 


The three test profiles consisted of two deep ocean 
profiles and one shallow water profile. The three profiles 
were the three test cases for the Acoustic Environmental 
Support Detatchment Workshop on Acoustic Propagation 
Modelling by Non-Ray-Tracing Techniques [8]. The profiles 
are listed in Tables 1 through 3. The results are plotted 
begining with Fig 16. 
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TABLE I 


AESD 1 PROFILE AND RUN TIMES 


PROFILE 
Depth Sound Speed Depth Sound Speed 
meters n/sec feet ft/sec 
0.0 1536.50 0.0 +: 5041.00 
152.40 to39. 24 500.00 5050.00 
406.390 1501.14 1333.00 4925.00 
1015.90 1471.88 3333.00 4829.00 
5587.89 1549.60 18333.00 5084.00 


Bottom Sound Speed: 1555.52 m/sec, 
Bottom Density: 1.9176 gm/cc 


5103.42 ft/sec 


Source Depth: 253.90 mn, goo U0 grt 
Receiver Depth: 863.50 a, 2833.00 £t 
Frequency: 25.00 Hz 


OMODE RUN TIMES 
Calculation of 44 Modes: 3.62 sec 


Transmission Loss to 120.0 nm 18.43 sec 


V7 





PROFILE 
Depth Sound Speed Depth Sound Speed 
meters n/sec feet ft/sec 

0.0 149. 12 O20 4984.00 
60.96 1520.34 200.00 4988.00 
91.44 7511.20 300.00 4958.00 
182.88 1507.24 600.00 4945.00 
487.68 100,08 1600.00 4934.00 
1219.20 1511.81 4000.00 4960.00 
1402.08 1508.76 4600.00 4950.00 
1828.80 1499.62 6000.00 4920.00 
2438.40 1504. 80 8000.00 4937.00 
3962.40 a2z7. 66 13000.00 5012.00 
5486.40 1554.48 18000.00 5100.00 


Bottom Sound Speed: 


TABLE II 


AESD 2 PROFILE AND RUN TIMES 


1560.39 m/sec, 


5119.40 ft/sec 


Bottom Density: 1.91/76 gm/cc 


800.00 ft 
3600.00 ft 


Source Depth: 243.84 a, 
Receiver Depth: 1097.28 a, 


Frequency: 50.00 Hz 


QMODE RUN TIMES 


81 Modes: 4.36 sec 


Transmission Loss to 300.0 nm 21.60 sec 


Calculation of 
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TABLE III 
AESD 3 PROFILE AND RUN TIMES 


PROFILE 
Depth Sound Speed Depth Sound Speed 
meters a/fsec feet ft/sec 
0.0 1524.00 0.0 5000.00 
15.. 24 1520.95 50.00 4990.00 
ie 12 onl? =) 150.00 4980.00 
Bottom Sound Speed: 1541.31 m/sec, 5056.80 ft/sec 


Bottom Density: 1.40 gm/cc 


Source Depths 6.10 om, DOS 00S Et 
Receiver Depth: 12.19 mn, 40.00 ft 
Frequency: 500.00 Hz 


OMODE RUN TIMES 
Calculation oft 5 Modes: 1.32 sec 


Transmission Loss to 25.0 nn 6.05 sec 
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In the preceding section approximate solutions to the 
homogeneous differential equation (2.10) were developed. [In 
this section a source term will be added to the differential 
equation, resulting in an inhomogeneous differential 
equation to be solved by a Green's function superposition of 
homogeneous solutions. 


A. SOLUTION OF THE INHOMOGENEOUS EQUATION 


The inhomogeneous version of eguation (2.10) for a 


source at depth z , and zero range, is 
Ss 
1 d_ (rd) + d26+ we =) 2a On rn 02 —Z Sal 
iE Sz | ar dz2 ce : Let Pe ! 
where 6 is the delta function. Applying the Hankel 


transform (Eas (2. 1%)) results in the inhomogeneous 


differential equation 
aa ck ier (z,K ) 26 ( ) (5.2) 
u + - u(Z = - Z=2e) - 
dz2 fe <: Ss 
The conditions imposed upon this equation are the same as 
(2.14) with the exception that gu is discontinuous at depth 
z 
w=Z First, the Green's function will be defined and then 


Ss 


the solution will be integrated using the inverse Hankel 


transforn. 
1. fhe Green's Function 
equation (5.2) can be solved by means of a Green's 


80 





function {1], defined by 


Me—mGtzaaey— —20 (Z)U0 (2 }) S W(U ,U ,z ) 
} 1 2 s 2s 


Oe Z Ss 2, (233) 
s 


and 


= G(z, ra u W(U ,U , 
u (z ae) a nae) v cae 3 z) 


S 
74 Ee eds 6 Sy 
= b 
The function U (2) is the solution to the homogeneous 
equation (2. 13) which satisfies the surface boundary 
condition. The function U (2) is the homogeneous’ solution 
which satisfies the lower boundary condition. Se ae is 
the Wronskian for the two solutions, given by 

ee = Sez) = Sie) 

where primes denote partiai differentiation with respect to 


depth. The separate solutions Ux(2) and oe 2) become 


Z 


~1/2 
U (z) =k cos( {k dz - 9), (5.4) 
1 Zz pez 1 
tall 
-1 atl 
U (z) =k cos([k az - 8), 
2 Z 


where Z eA I Ae 
t t 


and 8 and oe are determined by the surface and botton 
1 


boundary conditions. Designate 


Zz, 

GS fx qdz- 9, (5.5) 
Z 1 

Z. 

qo! 
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+1 
L = | k az- 9, 
Z Z 2 
z 
which gives an expression for the Wronskian, 
-1/2 


W= kK {d_ (LL )cosL sinL - 
dz 1 iz 1 


z 


d (L )cosL sinL 
Se | 2) 1 ai 


W = sin (L +L). Slesle 
1 2 


For those values of k for which Sin(L +L) equals zero, the 
r 


Green*s function is Singular. The characteristic equation 


for a zero Wronskian is 


tel 
L +L = k dz - 8 - 8 =nT 
1 2 7 Z 1 2 
pal 
N=O Fle 2, 343 6 (5.7) 
where n+1 124s the mode number. This is identical to 


homogeneous characteristic equation (2.30). The derivative 


of the Wronskian with respect to k is 
rc 


W' = dW = cos(L +L) d_(L +L ) 
dk™ eres we Oks le 2 
: -4 
= cos(L #L -k k dz - d (8 +@ 2 Sac 
cos { ; a 2 2 ( e aR 1 ae ( ) 
u 


Since W' is well defined at the zeros of the Wronskian, the 


poles of the Green's function are simple poles. 
2. Evaluation of the Integral 
TO evaluate the integral expression for (Z,C), 


6 (Z,r) = | seen ene I) CK, (5.9) 
0 s xr QO f ESSE 


it is standard practice to express J (k r) as a pair of 
QO oF 


Hankel functions (Ref. [2], [3]. [4]) 


a2 





C1) : €2) 
J (k = ; 
me a 4 ( H S19) He Sis 8 J 


Now becomes 
Es (1) 
(Z,r) = 1 G(Z,Z 3k) # (k rc) k dk 
220 Ss iE 0 r Er Ur 


i 2 
+ 1 fc ZZ 3K) ut (Kr) X dk. (5.10) 
as 0 s fr 0 rE cr UF 


Two properties of the Hankel functions are now used: 


1 
4M ). 


; QO for large argument in the first guadrant of the 


(2) 


complex plane, and 5 > Q for large argument in the fourth 


quadrant. The first integral is evaluated by integrating 


around the first quadrant, while the second integral is 
evaluated around the fourth quadrant (Fig 19). The poles of 
the integral lie on the real axis so the integrals in (5.10) 
are not well defined. However, the boundary conditions are 
met only by integrals whose path of integration in the 


complex plane passes below the poles. The real part of the 


horizontal wave number, k , represents the propagation of 
Is 


the wave functions with range. Any imaginary part of k 
represents an exvontential growth or decay of the wave 


function with range. For k to represent an outgoing wave 
r 


(taking into consideration the sign in Eq. (2.8)), any 


negative imaginary part represents waves which 
"exponentially grow" with range. Since this is inconsistent 


With the original restrictions on ¥, it is required that 
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Im(k ) > 0. This 1S physically appealing because it is 
rc 
eguivalent to having a small amount of attenuation. [In 


addition, there 1s a branch point at k equal nace where 
c 
c is the sound speed in the bottom. However, at long range 


the associated branch line integral has been found to make 
an insignificant contribution [5]. 
The first integral now may be represented by 
Lo 
(1) . (1) 
1 H (k cr) Gk dk + Tit Res{H (KP ECG ko} 4 (S611) 
2 Ty 0 iG im 65 0 c rc 


and the second by 
H (k rc) G k dk. (5.2112) 
r Eee 


Since 


(7) (2) 
LK =- =k 
H (ExeL} H (-LK> 2) 


and G(zZ,z ;k ) is an even function of k , the integrals in 
s or c 


equations (5.11) and (5.12) cancel. This leaves 


“ N 1 
(Z,C}) .= +71 f Res (H' ae EE) G2, Zack ) tes (54:13) 
1 Q r r Ss fC 


Since the poles of G(zZ,z ;k ) are simple and occur at the 
s fC 
zeros of the Wronskian, the derivative of the Wronskian with 
respect to k is used to find the values of the residues. 
rc 


Further, since at the zeros of the Wronskian, the two 


solutions et and ce are identical let 
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which yields 


° 1 
(Z,C0) = wr2ni 2 ra ik Eyes) U(Z yk | f Wt. (5.14) 
poles 0 rc See 


3. Normalization 
From the derivative of the Wronskian (Eq. (5.8)) a 
normalization factor is designated 
Z 
itil: 
[)u2q] [. tags aes ee 5.15 
u = z + + s e 
Zz 5 oat 1 2! 
Abu 
For long range, equation (5.14) can be written (using the 
asymptotic form of the Hankel function [6] 


: 1/2 az 
(2,0) = +2 (27) Ef (kK ©) exp{i(k ctr) } 
poles c E 5 


Oz 0 (2) / ilufil. (5.16) 


This use of the WKB approximation in the denominator of the 
Green's functiagn allows an estimate of the normalization 
factor, |jJu2i], to be made for the solution. This estimate 
allows calculation of the eigenfunctions at any particular 
depth, using only the WKB approximation. It should be noted 
that the skip distance of an equivalent cay, xX, is 
approximately equal to the normalization factor divided by 


the horizontal wave number, 
Zz 


cals oe 
X = k | k dz, 
C & 


teu 
> 


Seti] Ss 
c 
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Be. THE Q FUNCTIOW 


Pinally, a function Q (to be used in the next chapter) 


is defined 
Z 


Ca 
Q=if[ Ke a Ot er ye (5.17) 
ee u 1 
W Gu 


From the characteristic egquation (2.30) it is at once 
obvious that the mode eigenvalues correspond to integer 


values of Q. Differentiation of Q with respect to k 
ic 


provides a second form for the normalization factor, 


1 
piu2ti = k n | dQ / dk |. (5-18) 
c | rc 
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VI. PDHE MODEL 


In this chapter the computer model based upon the 
preceding theory is described. A later chapter will give a 
brief description of the actual computer program itself. 
The nodel performs five general steps in estimating 
propagation loss between two points: 

1. First the environmental information (sound speed 

profile, bottom type, etc.) is processed into ae forn 
Suitable for application of the preceding theory. This 


involves two major steps: the assumption of a sound speed 


profile consisting of layers withc varying linearly 
Within each layer; and a division of the profile into 
regions representing separate mode "families". A detailed 
description of the family grouping is given below. 


2. Next, the eigenvalues are determined by finding 


those values of the horizontal wave number k which yield 
rc 


integer values for Q (defined in Eq. (5.17)). The 


eigenvalues are found in two steps. First an initial 
estimation is made with a least squares polynomial 
estimator; then the eigenvalue is further refined with one 
to five iterative steps. 

3. After determining the eigenvalue for a particular 


mode, Kk , the value of the resulting eigenfunction 
Ge i 


U(z,k ) must be found at the source and receiver depths. 
gs 


In addition, a value of the eigenfunction normalization 
factor, |ju@|{|{, must be estimated (Eq. (5.15) and (5.18)). 


4. At this point, a set of functions (the norma 
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modes) has been defined which gives the depth-dependent 
behavior of the acoustic field. The next step is to 
represent the propagation of the modes, giving propagation 


loss as a function of range. This is done in two steps: 


a The modes are attenuated representing losses 
due to attenuation in the bottom, boundary roughness, 
and the empirical attenuation of sound in sea water. 

® The modes are summed including the coherence 
effects of the Hankel function of Eq. (5.16). 
is As the range between the source and receiver is 

increased, the source and receiver sound speed profiles 
may differ significantly, resulting in different sets of 
modal eigenfunctions. The transition of acoustic energy 
between differing environmental domains must now be 
treated. This is done employing the adiabatic assumption 


of mode invariance for slow horizontal changes. 


A. MODE FAMILIES 


To apply the WKB methods (especially the definition of 
Q) developed in the preceding chapters, the normal modes are 
divided into separate families based upon their position in 
the sound speed profile. The sound speed profile is divided 
into various domains by depth and phase speed (or horizontal 
wave number), each domain representing a mode family with 


particular characteristics. 
1. Definition of Families 


> > aes oe sis ee ei ae ee 


If there exists more than one possible combination 


of upper and lower turning points for a particular k (that 
r 
is, there is more than one acoustic channel), the definition 


of Q becomes ambiguous. The ambiguity in Q , and the 


resulting ambiguity of modal eigenvalues, is resolved by 
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grouping Q (hence the modes) into separate families 
depending upon the position of the turning points in the 
sound speed profile. 


In the bilinear profile of Fage20, for k 
iE 


corresponding toc ; the definition of Q (Eq. (5.17)) may 
Pa 


be be applied for turning points either between 2 and 7 or 
oO 


between Zr and ee The separate applications represent 
Separate interpretations of Q and separate families of 


modes. This profile may be divided into three families. 


The upper family (I) 1s limited by o = <¢ S$ a and 
P 


f= Z eZ $< z. Modes in this family correspond to waves 
O te )}6hCU 1 


which are refracted at the lower turning point (above a 


and reflected at the surface. The lower family (II) is 


lA 


limited Dyers Cc S CC and Zz. SZ. ,Z 
3 Pp 2 1 cu td 


this family represent waves which are refracted at the upper 


Ze Modes in 
b 


turning point (below 2m) and reflected from the bottom at 
depth ae For phase speeds above = the Q function and 
modes are unambiguously defined, and a third family (III) is 
delineated byc Sc S$ c and2z2 ‘SZ ,zZ < Boe Modes in 


2 p b O ae 
this family represent waves which are reflected from both 


the surface and botton. 

For this simple profile three families are 
delineated by both the local maxima and minima of the sound 
speed profile, and the surface and bottom boundaries. In a 
Similar manner, a sound speed profile of any complexity can 
be divided into separate families: for every sound speed 
Maximum two separate families exist with turning points 


above or below the depth of the sound speed maximun. 
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FIGURE 20 - Bilinear Profile Families 
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2. Family Types 


The sound profile is divided into separate families 

BOL : multiple channels, purely refracted waves, bottom 

reflected waves and surface reflected waves. Illustrations 

of the family classifications are given in Fig 21 and 

Fig 22. The purely refracted modes are divided into four 

types: 

a Type 1- refracted - refracted (RR). A mode which is 
refracted at both the upper and lower turning points. 

e Type 2 ~ RR with upper channel. A mode which is 

refracted at both the upper and lower turning points. [In 


addition, at the minimum k for this family, there is at 
rc 
least one family above it in depth. 


a Type 3 = RR with lower channel. A 
refracted-refracted mode for which there exists at least 
one family below it. 

e Type 4 = RR with upper and lower channels. A 
refracted - refracted mode for which there exist families 
both above and below it. 

In a Similar manner the modes which are refracted at the 
lower turning point and reflected at the surface (refracted 
- surface reflected or RS&) are divided into two types: 

a Type 5 - RSR with no other channels, and 

e Type 6 — RSR with lower families. 

The modes which are refracted at the upper turning points 
and reflected at the bottom (RBR) are also divided into two 
types: 

ea Type 7 - RBR with no other channels, and 

a Type 8 —- RBR with at least one upper channel. 

Finally, we have the last family, 
e Type 9 - SRBR, modes which are reflected by both the 


bottom and surface. 
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PIGURE 21 - Family Classifications I 
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FIGURE 22 - Family Classifications II 
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ae Single Channel 
Within a single channel profile, as _ changes 
from one family to the next, Q is a continuous function of 
oF since i 0G and the integral of the horizontal wave 


number as functions of k are all continuous. Thus, the 
c 


upper limit of Q(k ) for one family is equal to the lower 
limit of Q(k ) for the next. Although this result may seen 
obvious, it is not generally observed with a multiple 
channel profile. | 

b. Double Channel 


The WKB methodology used in this model treats 
the modes in the upper and lower channels of a multiple 
channel profile as separate and independent phenomena. As 
long as the barrier between the two channels is of 
sufficient width this 1s an accurate assumption; the sound 
fields in the two channels remain almost independent of each 
Other. However, as the barrier decreases in width, the two 
sound fields no longer remain independent; the modes in one 
Channel become increasingly coupled to the acoustic field in 
the other channel. When this occurs the WKB methodology 
begins to break down. 

Consider two families separated by a thin 
barrier (Fig 23). For modes with phase speed close to the 
edge of the barrier, the integral of the imaginary vertical 
wave number across the barrier (L in Eq. (3.15)) approaches 
zero. As this occurs, the value of the lower turning point 
phase angle in the upper channel is perturbed by increasing 
amounts away from the base value of r/4 (Eg. 3.16)). The 
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size and direction of this perturbation depends upon the 
oscillatory solution in the lower channel. 

If there is a large number of closely spaced 
modes in the lower channel, the WKB solution at the botton 


of the barrier oscillates rapidly as a function of k , which 
c 


in turn may produce rapid oscillations in the value of the 
lower phase angle in the upper channel. If, in addition, 


the integral of c in the upper channel varies slowly with 
at Q(k ) changes in Character from monotonic ZO 
oscillatory. If this happens, it is possible for more than 
one value of k to satisfy the WKB characteristic equation 


Cc 


in the upper channel, and there occurs the possibility of 


multiple solutions for a particular mode. It should be 
noted that for profiles based upon actual oceanographic 
conditions, the integral of the vertical wave number across 
the channel usually increases rapidly enough to ensure that 


O(k ) is monotonic. 
rc 


A second difficulty occurs at the crossover 
between the two multiple channels and the main channel 
meng 2 ). The Nth mode in the main channel represents the 
Nth possible solution (in order of increasing phase speed), 
requiring that exactly N-1 distinct solutions exist at lower 
phase speeds. The maximum number of modes in each of the 
upper and lower channels is the largest integer less than 


the respective maximum Q(k ). The first mode in the main 
E 
Channel is the integer larger than its minimum Q. Assume N 


is the first mode in the main channel; for N-1 solutions to 
exist at lower phase speeds, the sum of the integer parts of 
the upper and lower channel maximum Q's must be N-1. FOr 


the upper channel, the maximua Q is 


oF 


Q = 
ite 


ai 


ei 
( | Kode Oe GS oT de 
Z Z u 1 
tu 
for the lower channel, 
Z 
wall 
Q = 1 kK dzZ=-6¢6 ~=- 89 +1 
II rt | Zz 2 i Je 
2 


and the minimum Q of the main channel is 


eel 
Kode = Oe Oe jc 
Z u 1 


teu 


The above expressions for Q , Q and Q show that the sum 
I If Lit 


Q = if 
avec = 


of maximum modes in the upper and lower channels may be 
equal to N~2, N-1 or Nw Thus it 1s possible for the WKB 
characteristic equation to have no_ solution or multiple 
solutions for a particular mode near the crossover between 
multiple channel and main channel families. This phenomenon 
is a mathematical consequence of the fact that near the edge 


of the barrier, both k and g are sufficiently small to 
Z ° 


invalidate Eg. (3.17). 


This difficulty could be overcome in a number of 
Ways, which‘include changing the nature of the barrier or 
using a method other than the WKB method in the region near 
the barrier edge. If the edge of the barrier is 
appropriately smoothed so that the derivative of the 
vertical wave number is proportional to the square of the 
vertical wave number, the WKB validity criterion will remain 
satisfied. For this particular model this scheme would be 


undesirable, as it would require the use of a profile other 


ee 
than the c linear profile. The c profile is used 


primarily for the rapidity with which the vertical wave 
number may be integrated. Since integration of the vertical 
Wave number is central to many Operations in the 4KB method, 
this integration time is extremely important to the progran 
execution time. Another method to overcome this difficulty 
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would be the alteration of the barrier connection formulae 
(Eq. (3.15) ) for very small values of the integral L. A 
third method would be to modify the integral of the vertical 
wave number, so that a correction term is applied at depths 
which fail to meet the WKB validity criteria. Again, this 
method adds to the execution time required for the 
integration loop. In the program described in the next 
chapter the difficulty of a duplicate eigenvalue is overcome 
by eliminating the duplicate eigenvalue from the main 
channel. The difficuity of a missing eigenvalue is overcome 
by assuming the missing eigenvalue to be at the horizontal 


wave number corresponding to the edge of the barrier. 


Be ELGENVALUE ESTIMATION 

Following definition of the mode families, the next step 
is the determination of the mode eigenvalues. 

1. fhe Eigenvalues Search 


The search for eigenvalues is done by families, 
determining those values of the horizontal wave number which 
most nearly yield integer values of Q. The eigenvalues are 


found in three steps: the formation of an initial 


polynomial estimator; the finding of two values of k which 
c 


bracket k ; then performing successive iterations by the 
c,u 


method of regula falsi. For each family, sample values of Q 


as a function of phase speed are found (by integration of 


the vertical wave number and determination of oF and @). 
u 
Through these points a fifth-order polynomial in Q(K ). 


p(Q), is fitted by the method of least squares. This 
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polynomial is used to find the first estimate of the 
horizontal wave number. For long ranges and when 
considering incoherent transmission loss, this estimate is 
sufficiently accurate, and the subsequent iterations may be 
bypassed. Using the derivative of the polynomial estimator, 


p'(n), a second estimate for k is then sought which, in 
Sei 


conjunction with the first, brackets the eigenvalue. From 


these two estimates which bracket the eigenvalue, the 
eigenvalue estimate is improved by successive iterations, up 


to four times (the method of regula falsi [{1]). 


2. Selection of the Estimator 





A number of techniques were considered as estimators 


for k as a function of Q, including a cubic spline and the 
rc 


least squares polynomial (which was used). 


a. Cubic Spline 


The cubic spline fits a piecewise cubic 
polynomial through a given set of data, resulting in a 


function continuous through the second derivative. When 
Q(k ) was monotonic and reasonably smooth the spline gave 
c 


excellent results (in fact, it was exact for certain 
analytic profiles). In the case of the multiple channel 


profiles, QO(k ) may become irregular and slightly 
c 
oscillatory in character. .In this case an unfortunate 


combination of sample points wiil generate a spline which 


gives estimates of k in the wrong order (Fig 24). The 
r 
cubic spline was deemed too sensitive to the minor phase 


effects of the muitiple channei modes, and a smoother 


substitute was sought. 
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FIGURE 24 - I1l-Conditioned Cubic Spline 
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Ds Least Squares Polynomial 


The problems with the cubic spline ied to an 
estimator which would be insensitive to any small _ scale 
oscillation in the data, yet which reliably estimated the 


general trend of Q({k). A least squares polynomial 
LE 
algorithm using orthogonal polynomials from the text by 


“Conte and deBoor [2] was selected. To decrease the 
difficulties near the edges of multiple channel barriers, a 
weighting was used which de-emphasized the minimum and 
maximum values of Q in the families. The polynonial was 
found to be of sufficient accuracy to provide, in most 
cases, an estimate accurate to within one part in twenty of 
the intermode spacing. In addition, the polynomial provides 
a rapid estimate of the derivative of Q(k ) with respect to 


k 3; this result is useful in later iteration steps and in 
rc 


normalizing the eigenfunctions. 


C. CALCULATION OF THE EIGENFUNCTION 


After the eigenvalue k is determined, the next step 
c,n 
is to calculate and normalize the eigenfunction JU (z) at the 
n 
source and receiver depths, Z and Zz. To determine the 
Ss iE 


eigenfunction, one of the three solutions developed in the 


theory chapters (the WKB and Airy solutions), or a_ function 


extrapolated from the known solutions is used. 


1. Assumed Solutions 





At the receiver and source depths, the Airy and WKB 


validity criteria, F and P , are evaluated. If a WKB or 
Ww a 
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Airy solution is a good approximation, the eigenfunction is 


determined in two steps. 


a The coefficients for the WKB solutions, a,b,c, 
O O O 


and d are determined for each ensonified and shadow zone 
O 
region, using the values of the upper and lower turning 
point phase angles 8 and oe 
u 
wi If the WKB solution is a good approximation, the 
integral of the vertical wave number from the upper or 
lower turning point to the source or receiver depth is 
determined, and the eigenfunction is evaluated using 
Peet 2.22) or (2.26). 
e If the Airy solution is a good approximation, the 
Airy coefficients, A and B, are determined from the 
O O 


connection formulae (Eq. (3.13));3 the variable x is 


calculated; and the Airy solution is evaluated with 


Eg, (3.5). 
2. Spliced Solution 


If neither the Airy or WKB solution is a_e good 
approximation at the desired depth, a finite difference 
technique is used to extrapolate the eigenfunction from the 
two accurate WKB oor Airy solutions above and below the 
desired depth. A search is made from the desired depth to 


find the nearest depths above and below, - and ae for 
which either the Airy or WKB forms are good approximations 
Mig 25). At these depths oe) and otze) and their 
derivatives with respect to depth, ue and ea? are 
calculated. The depth interval between Es and oe is divided 


into (up to 100) equal increments, z_, and a solution is 
i 
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iterated from aan to ze using a finite difference equation. 
The extrapolated function is normalized to eG then the 
procedure is repeated from = to a The final 
eigenfunction is a weighted average of the two extrapolated 
functions. 


3. EHigenfunction Normalization 


Based on the preceding theory, two methods of 
normalizing the eigenfunctions are available: evaluation of 
the integral of Eq. (5.15); or estimating dQ/dk from the 

rc 


polynomial estimator to use in Eq (5.18). In the case of an 


eigenfunction for which there are no upper or lower 
channels, Eq. (5.15) is evaluated. When an upper or lower 
Channel exist at the same phase speed as the eigenfunction, 


the evaluation of the terms for 98 ' and oo becomes 
u 
difficult to evaluate across the barrier. In this case the 


normalization factor is evaluated with the derivative of the 


eigenvalue estimation polynomial, and Eq (5.18) becomes 


-1 
fiuei] = Jk dQ/fdk {| = c k (p'(Q))  .« 
Ee i Dp, 0) 5, a 


In the sections that follow, u (z) will represent the 
n 


normalized eigenfunction, 


72 
u (Zz) = U (z) / (1luett) . 
n n 


D. TRANSMISSION LOSS 


The object of this model is the calculation of 
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transmission loss between some fixed location and a second 
location which varies in range from the first. Inspection 
of the equations used herein will reveal complete symmetry 
with respect to the source and receiver locations. Thus the 
labelling of source and receiver locations is completely 
arbitrary, irrespective of which location may in fact be the 
site of generation of acoustic energy (the principle of 
acoustic reciprocity). Accordingly, in what follows the 
source location refers to the location which is fixed in 
Space, and the receiver refers to the location which varies 


in range. 


1. Coheren 


ct 
{tJ 
in 


ansgission Loss 





The transmission loss is defined by 
2 2 
TL = 10 log {P /P {, (6. 3) 
10 O ic 


where P andP are the acoustic pressures at unit distance 
oO ie 


and range r, respectively. Since the density of seawater 1s 


assumed constant, and the velocity potential has been 
normalized to one at unit distance, Eq. (6.3) reduces to 


2 
TL = -10 Boor, ld |- (6.4) 


At this point the effects of attenuation of the individual 
modes due to absorption of sound in seawater, botton 
absorption, and scattering at the bottom and surface must be 
included. At frequencies below 1000 Hz, the absorption 


coefficient is given by the empirical equation [3] 


2 
a= 125 £ , (6.5) 


where a is in decibels per nautical nile, and f is in kHz. 
This tern, the effects of absorption of sound in the botton, 


and effects of bottom and surface scattering will be 
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combined into a single term, D, defined in a later section. 
n 


With this attenuation term, a mode pressure term is defined 
~1/2 
P = (k r) u (Z ) wu (Zz ) exp(-D r); (6.6) 
n z n s n ox n 
Eq. (6.4) may be re-written 


2 
TL = -10 log. { % P cos(k r)j + 
10 nh a c,n 


f 


2 
{ = P sin (k Eels (6.7) 
nh on en 


¢ 


2. Averaged Transmission Loss 


The transmisson loss of Eq. (6.7) is the sum of 
acoustic energy in each mode plus the effects of 
interference between modes. The resulting transmission loss 
1s characterized by fluctuations with range, the rapidity of 


which depends upon the differences between the eigenvalues 


AE = 2 f Ak e 
n,m r,,m 
where Ar 1s the range of the fluctuation and 
n,@ 
Ak = (k ee he (6.9) 
r,n, c,h c,m 
For adjacent modes (n=mt1), Ak depends upon the gross 
c,n,m 


features of the sound speed profile. For widely separated 
modes the difference between eigenvalues becomes dependent 
upon the smaller details of the sound speed profile. Thus 
due to irregularities in the medium and the inaccuracies of 


the model, the higher values of Ak are estimated with 
r,n,m 


less confidence, and the more rapid fluctuations (especially 


at longer ranges) are not considered an exact representation 


of the sound field. These rapid fluctuations can be 
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smoothed by faltering the mode interference wave numbers 
corresponding to the shorter wavelength fluctuations 
Ak Dee e/ oh 
C,0l,m Ss 
A cunning average, which is a window in range space, 
represents a filter in wave number space. Thus a running 
average of the acoustic pressure squared term, optionally 
weighted with a Gaussian-.filter, is used to give an averaged 
transmission loss figure. The Gaussian filter was selected 
since the Gaussian weighting function in range transforms to 
a Gaussian weighting function in wave number space, allowing 
a smooth roll-off from long wavelength to short wavelength 
fluctuations. This physically corresponds to the 
elimination of the finer details of the sound speed profile 


aS contributors to the calculated acoustic field. 


3. Incoherent Transmission Loss 


== Gat ae Sa ce Se ec eb Se oe aD 


At long range ("long" depends upon the acoustic 
frequency and the inhomogeneities in the ocean), the 
coherent interference effects lose significance and the 
incoherent acoustic pressure (representing an RMS sum of the 
energy in the modes) becomes the most useful estimate for 
the acoustic field. Accordingly, incoherent transmission 


loss is defined 


TLI = -10 log EoD (6.10) 
10 om a 


Ee MULTIPLE PROFILE PROPAGATION 
In the preceding sections no horizontal variation in 


both sound speed and bottom depth has been assumed. Within 


one locale (at sub-kilohertz frequencies) this assumption 
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can be fairly accurate. However, as sound propagation over 
greater ranges iS considered, markedly differing oceanic 
regimes are encountered, and this assumption becomes 
untenable. Thus the long range sound propagation model must 
include the effects of changes in the sound speed structure 


With range. 





1. Adiabatic Assumption 


TO incorporate horizontal change the so-called 
adiabatic assumption is invoked, namely: the energy in any 
particular mode stays within that node. The assumptions 
involved may be more explicitly stated. TO paraphrase 
Tolstoy and Clay [4]: 


»s The eigenfunctions and eigenvalues are dependent upon 
the local stratification. 

es The acoustic medium varies slowly with range. 

»e There is no cross coupling between different modes in 
the transition between profiles. 

a There is no back reflection of acoustic energy at the 


transition between profiles. 


Milder has shown that the accuracy of the adiabatic 


assumption decreases with increasing frequency and 


4 (a) (b) 
horizontal gradients of k [5]. iu and u are 
m n 


adjacent modes (nm+i), Milder shows the adiabatic assumption 


is a good approximation if 


Z 
1 (a) 2 (b) 2 2 
(kX) X { fu d <«k ) u dz} / k << 4m, (Oe) 
wv dr n 
tu 
where X is the skip distance for an equivalent fray. 


Milder's equations may be extended to estimate an order of 


magnitude transition range required for an adiabatic change 
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between two profiles. If the magnitude of the integral in 
Eq. (6.11) 1S eStimated by 
Zz 
aaa) 2 (b) 
u _(k ) dz 
7 m c 
tu 
» 2 e 2 
= d(k dz = 1 Ak dz 6.12 
ap efx)! fr = | ‘ ( ) 
and an expression for Ar, the distance required between 
the two profiles for an adiabatic change to occur, is found 
2 2 Z 
NED Dek [Ak Jdz Y/Y (4/2 7 ). (6.13) 
This transition range is calculated by the program and 
printed out as a guide to the applicability of the adiabatic 
assumption. 
3. Iransmission Loss 


Since the eigenfunctions and eigenvalues are assumed 


to be characteristic 


of the local profile, the source and 


receiver mode excitation is expressed in terms of the 
, (S) (rc) 
eigenfunctions at the source and receiver, u and u 
n 
respectively. Defining 
K (c) fr dr, (6.14) 
n ie, 0 
0 
and 
oo rc 
D (r) -f De dry, (6.15) 
n n 
0 
the adiabatic mode pressure term is derived, 
= 1/72 (S) (r) = 
P = 2[2 1n/K (r) ] u u exp(-D (r) ]. (6.16) 
D m m m m 
This gives the transmission loss for the multiple profile 
case: 





= 2 
TL = -10 log el Pecos (kK (r)) ] ~ + 
1 n a n 


- 2 
ire Sin tn e(r)) } . (6.17) 
bh a Mm 


4. Mode Matching 


Multiple profile propagation using the adiabatic 
assumption requires the matching of modes between profiles. 
In previous works this has been done by numbering the modes 
in order of increasing phase speed. This is always correct 
for the single channel profile, in which the mode order is 
always invariant. However, with multiple channel profiles 
this procedure may yield physically unreasonable results. 

Consider two sound speed profiles with two near 
Symmetric channels, and six modes (Fig 26). Since the upper 
channel is slightly "deeper", the modes in this channel 
occur at a lower phase speed than the corresponding modes in 
the lower channel. Thus, in terms of overall order, the 
odd-numbered modes represent acoustic energy trapped in the 
upper channel and the even numbered modes represent acoustic 
energy trapped in the lower channel. Slowly alter the 
profile so that the lower channel becomes slightly "deeper" 
than the upper; and the modes in the lower channel precede 
the modes in the upper channel. Matching modes by overall 
order would, in this case, couple all the modes in the upper 
Channel with modes in the lower channel, and vice versa -a 
physically implausible event. In matching between profiles, 
not only mode order, but mode families must also be taken 
into account. 

Assume again the profile changes slowly, this time 
so that both upper and lower channeis become "shallower" 
1g 27). Two modes appear in each of the upper and lower 
sub-channels, and two modes appear in the main channel. Tt 
is not clear which of the new nodes in the main channel 


corresponds to the missing last modes from the upper channel 





SOPpOW YATRA STouUeYD OTAZONNAS ARON OAL — 97 ANNOTA 


Si9ona 19nd 19nd 
,SaMOTIVHS.., YAMOT ,Y3ad3930, J Uadd/N ,Y3d33G,, 


\ 
is it 
: rp ere | : . 


peeds punos 





Yidep | 


112 





SUOTITSURIZ TauUeyD Sssory - ¢z TANNTA 


IVAIOYd WMA 


Ng W4dOUd ONINNIDAE 


waeceryrcre-- 






Yidep 








SI 140d 
ILVIGSNYILN 





113 








and lower channel respectively. In fact it is not clear 
that the transition from one family type to another is an 
adiabatic process [6]. Clearly some consideration must be 
made for the uncertainties of mode transition between 


channels. The matching of modes proceeds in two steps, 
a. Same Channel 


Modes within the same channel of the sound speed 
profile (main, upper, middle, lower or surface ducts) are 
matched first. This is done by searching the modes for both 
profiles in order of increasing phase speed (lowest modes 
first), matching nodes with the same mode number and in the 
same channel. In most cases, this provides a amatch for 
nearly all the modes. If one profile has fewer modes than 
the other, the highest modes are truncated, physically 
corresponding to leakage into the bottom. At this point 
almost all (all in the case of the single channel profile) 


of the modes have been matched. 
b. Cross~Channel Transitions 


The second step is the matching of modes which 
undergo a transfer between different channels. The results 
of such a transfer of acoustic energy depend upon the manner 
in which the sound speed profile changes with range. 
Consider the transfer of modes from both an upper and lower 
Channel into a main channel (Fig 27). Such a transfer could 
take place in two hypothetical steps, the order of which 
determines the outcone. 

First, assume that the lower channel becomes 
"shallower" and the last mode in the bottom channel is 
adiabatically transferred to the main channel. 
Subsequently, the upper channel becomes shallower, and the 
last mode in the upper channel is transferred to the main 
channel. In this case, the lowest numbered mode in the main 
channel corresponds to the last upper channel mode and the 


second main channel mode corresponds to the last lower 


channel mode. If the order in which the channels are 
perturbed is reversed, the order of mode transfer and the 
mode correspondence is reversed. 

Since the only information available is the two 
sound speed profiles, this situation results in an 
uncertainty in the matching of modes between profiles. 
(This uncertainty is analogous to the uncertainty pointed 
out by Smith [7] for adiabatic ray coupling between 
profiles.) A Similar uncertainty occurs when more than one 
mode 1s transferred from an upper or lower sub-channel to 
another sub-channel. Given exact information concerning the 
sound speed profile changes with range, the uncertainty can 
be resolved by computing the cross-coupling between modes. 
Milder has given the equations for the calculation of cross 
coupling between modes. However, these eguations involve 
the inner product integral of BGs ccs it), and the 
evaluation of this integral is both ill-suited to the WKB 
solutions used in this model, and too costly in terms of 
program size and execution time. Accordingly, a simpler 
formalism, similar to that used by Smith for ray-~tracing, is 
required to extend the adiabatic assumption to the case of 
multiple channel profiles. Thus the following rules are 
asserted to govern the transfer of nodes between channels: 


ls The total energy in the modes is conserved 


Zu. (Z ; 6.18 
an Beatie Z ( 


Zu ae 
hn s |profile 1 
ne If no uncertainty exists, that is if a transfer 
between a main channel and only one sub-channel 18S 
involved, the transfer will be assumed to be adiabatic 


and the modes matched by increasing phase speed. 


Bie If an uncertainty does exist, the energy of the modes 





involved is combined and then is divided evenly among 


the modes. 
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In physical terms these rules state that the 
mode transfers between channels which involve no uncertainty 
are considered adiabatic. When an uncertainty does exist, 
the exact nature of the mode transition between profiles is 
undetermined, and the mode energy is spread among all the 
modes undergoing a Similar transition. Although imposed as 
a result of the lack of information about the exact process 


taking place, the physical consequence is not unrealistic. 


F. BOTTOM MODELS 


In the preceding theory the boundary condition at the 
bottom was represented by the parameter we. In this section 
the various bottom models are developed which give 
expressions for both vu and the mode attenuation due to 


bottom absorption and surface and bottom roughness. 


f#. Fuli Fluid Botton 


The most commonly-used bottom model (1n normal node 
calculations) is a fluid half-space of density At the ocean 
bottom the condition of continiuity o of acoustic pressure 
and vertical particle velocity (Eq. (2.14) cand d) are 


imposed. The solution, Hac to the vertical Helmholtz 
equation in the bottom may be cepresented Dy the WKB 


exponential solution (Eq. (2.26)). With application of the 
continuity conditions, and reguiring that the solution decay 


at great depth, 


= = - ; 6.19 
Bee eo) Utz) ooo meee.) ( ) 


where 


Don 2 172 
B= [k -w fc] . (6.20) 
E 





Thus the full fluid boundary condition is written 


u = - p/o Bw (G_217) 
Oo b b 


The derivative of wu with respect to the horizontal wave 


number 1S required for the evaluation of cae in Eq. (4.19), 
tee Ok) 7) fp 8). (6.22) 
Oo.L b »b 


Given an attenuation coefficient, one for the botton 


Material such that 
Ge Cae 
b iG 


the mode attenuation due to the effects of an absorptive 
bottom may be estimated by [8] 
Zz, 0 
2 o2 2 
a =a 0 Jo 9GZ 7 aleo UrdZeet sp U dzj. (6.23) 
m b bY »b O bs 6b 
D O D 
The denominator of the above expression is related to the 


normalization factor, |ju2i], by 


Z, 


b2 7 
Jiju2fj = 2 co fo Zs o fo, dz j. (6.24) 
0 


or 


Substituting and integrating the solution in the bottom, the 


mode attenuation coefficient may be written, 


Z 
7 oF Ce (2) , aoe A LU, 


2. Rigid Bottoa 


The rigid or impenetrable bottom requires that _ the 


vertical particle velocity vanish at the botton, or 


p= 0, (6.26) 


ty 





and 
tea O (6.27) 


Since the acoustic wave does not penetrate the bottom, it 
is unaffected by the bottom material, and the rigid bottom 


is considered lossless, 


0. (6.28) 


RQ 
5 
it 


3. The Modified Fluid and Rigid Bottoms 


ae Ss ee a= — —=_ aan oo SSeS == 


In the modified fluid and rigid bottoms, the 
effective bottom losses are estimated by input bottom loss 
VS grazing angle data. The boundary condition at the 
bottona, uu, is calculated in the same manner as the full 
fluid and rigid cases. Given the bottom reflection loss per 
bounce, BL (in decibels), as a function of grazing angle, 
the attenuation coefficient per unit range can be written, 


fe 2.502905) BL /(20X), (6.29) 
n 


where X is the inter-bounce or skip distance. 
4. Bottom and Surface Roughness 


Terms are included in the mode attenuation for 
scattering due to a rough surface and bottom. From Tolstoy 
and Clay {9] an expression for mode attenuation due to an 
irregular boundary is used for the surface 

Zee 


2 2 
a = 2k S / X, for 2k. -S << 4 (6.30) 
s Z Z 


where S is the surface roughness. Similarly, for the botton 


Ze 2 2 
a = 2k B / X, EOQrmak <<. t, (6.31) 
b Z Z 


where B is the bottom roughness. The total mode 


attenuation is the sum of the attenuations due to bottom and 


118 





surface roughness, bottom absorption, and the absorption of 


sound in seawater, 


Denson to.) teem + Co (6532) 
m m b Ss 
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VII. THE PROGRAM 


ae ee eee SS 


This chapter contains a description of the computer 
progran, QMODE, which implements the model. The 
description, which includes a summary of the main progran 
and subprograms, is intended primarily for the user. QMODE 
is written in IBM FORTRAN IV(G), and was developed on the 
Naval Postgraduate School IBM 360/67 computer. The total 
progran Size (including object code, storage, and I0 
buffers) is 150K eight-bit bytes (170K for offline 
plotting). 


A. MAIN PROGRAM 


The main program is executed by cuns and profiles. A 
new run is started for each new source location; successive 
profiles within a run represent a succession of sound 
conditions with range. 

The input data are read uSing the NAMELIST feature of 
IBM FORTRAN, which allows the setting of a large number of 
default values to be optionally overridden by the user. The 
profile is listed in a vector PROFIL which can be tabulated 
in one of two basic forms: 

x Depth vs sound speed pairs (listed in either meters 
or feet). 

a Triplets of depth (in either meters or feet), 
temperature (in either degrees Farenheit or Celsius), and 
Salinity. 

The profile is then transformed into values of depth vs. 
sound speed in meters by using the appropriate conversion 
factors and Wilson's eguations for the speed of sound in 


seawater [1]. If the last sound speed vrofile point is 
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Shallower than a designated depth (DEPTH), the profile is 


extended to the bottom in an isohaline, isothermal gradient 
—1 2 
of 0.017 sec . Next the values of k and k for each depth 


point, and the values of y (Eq. (2.2)) for the layers 
between depth points are calculated. 

Calls to three subroutines complete the preparations for 
calculating the eigenvalues and eigenfunctions. A call to 
BOTSET sets the parameters reguired for calculation of the 
bottom boundary condition, u , and the mode attenuation 
coefficients. The subroutine BCSET sets parameters for the 


calculation of the upper and lower turning point phase 


angles, 94 and oa Next a maximum of ten families are 
u 


found, and their characteristics listed by the subroutine 


FINPAM. 
At this point the calculation of the modes begins in two 


nested loops: a family loop anda mode loop. For each 


family a sample set of c and the corresponding values of 
P 
Q(k ) are found by the the subroutine QFIND. A fifth order 
ie 


least squares polynomial, p(Q), is then formed through the 


sample points by the subroutine ORTPOL. 

Within the mode loop, first the eigenvalues are found; 
then the eigenfunctions and mode parameters calculated. fhe 
tirst estimate of a particular eigenvalue, k » is found by 

coral 
calling the polynomial estimator ORTGET, which returns an 
estimate of c for a given Q, 
P,n 


cl) = p(Q=n), 
Pp, 


a 





giving the first estimate of the horizontal wave number, 


k€2d = w/ cl), 
ea p,D 


for which a corresponding value of Q is found, 


QC1) = Q(kC1)),. 
aon 


The second estimate of k is then sought which will 
c,n 


"overshoot" the eigenvalue (to have two estimates which 


bracket the eigenvalue). The subroutine ORTGET also 
provides an estimate of the derivative of Q with respect to 


Kk, 
c 


-1 
C70 Pace (D010) ae. 
r Pp 


then, using a relaxation factor F to "overshoot" the 
c 


eigenvalue, the second guess is found, 


kC2) = kC1) + FP (QC1I=n) dQ/dk . 
Tat ees c c 


F is initially set at /2 and expanded by powers of J/2 until 
iE 
a value for k€2) is found for which, 
r,D 


(Q(k C22) =n) *(Q(kCt2)-n) < 0. 
,o aon 
Once two bracketing values are found, the iteration 
continues using the method of regula falsi until one of two 
conditions is met: a value of Q within 0.0005 of the 
ijateger nis found; or a total of five iterations have been 


made. Once the eigenvalue k has been found, a number of 


mode parameters are then calculated and stored for later 
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use. These include: k (KRVEC), mode order (MODENO, 
cen 


MORDR2), the mode number (MODEN), relative position in the 
sound profile (IBRN2), the skip distance (SKPDIS), and 
O(K ie (QN). The subroutine EIGFUN caiculates the value of 

’ 
the eigenfunctions at the source and receiver depths. The 
eigenfunctions are then normalized, completing the mode and 
family loops. 

If the current profile is the continuation of a multiple 
profile run, the current node set is matched with the source 
location modes by the subroutine MATCH. The subroutine 
PLOSS calculates the transmission loss from the range RSTART 
to RSTOP in increments of DELR. If the next profile is to 
be a continuation of the same run, two subroutines are 
calied to prepare for the continuation. The subroutine 
PLSET attenuates the source eigenfunctions and keeps a 


running sum of K (rc) (KRERV); subroutine SWAP keeps track of 
n 
the matching parameters. If, On the other hand, the next 


profile is the beginning of a new run, the default 
parameters are re-set. fhe input data are then read for the 


next profile and the main program sequence begins again. 


Be. SUBPROGRAMS 


1. BOTSET 





The subroutine BOTSET selects the bottom type (IBT), 
sets parameters for the calculation of the bottom boundary 
Parameter, uu, and the mode attenuation. The functions EMU 


and EMUP compute u and ue 6 as)6h6a: C6oEunction of k. The 
E 
function DECAY computes the mode attenuation as a function 


of the bottom type and characteristics, bottom roughness 


123 





(BOTRUF) and surface roughness (WAVEHT). If the modified 
fluid or the modified rigid bottom models are selected, 
DECAY linearly interpolates the bottom loss values from 


input data of bottom loss vs grazing angle. 


2. BCSET 


The subroutine BCSET divides the bottom and surface 
boundary conditions into phase speeds based upon the values 
of the Airy and WKB validity factors, F and F, versus the 

a W 


parameters KALIM and KWLIM. For the surface boundary the 


WKB forms are used at phase speeds below CP3 and above CP1; 
and the Airy forms are used at phase speeds between CP1 and 
CP2. For the bottom boundary the WKB forms are used at 
phase speeds below CP?7 and above CP4; the Airy forms are 
used at phase speeds between CP5 and CP6. FOE sthes "gaps 
between these phase speeds (two at the bottom boundary, one 
at the surface), third-order polynonials are formed to 


provide continuous functions for 8 and 8 . The subroutines 
u 


SBC and BBC calculate the values of 8 and oo as functions 
u 


oe K In order to distinguish between multi-channel 
rs 


families, a minimum Or maximum family depth is reguired for 


SBC and BBC respectively. When the parameter JP is set to 
-1 the subroutines also calculate values for 8 ! and a 
‘ 7 


3. FINFAM 


The subroutine FINFAM defines each mode family by 
checking for family boundaries at the surface and botton 
sound speeds, and each relative sound speed maximun. Two 
arrays are used to store the parameters for each family: 


maximum k , minimum k , minimum Q, and maximum Q (in _ FAM) 3 
E c 


and minimum depth grid, maximum depth grid, family type and 
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relative position in the sound speed profile (in IFAM). 
4. ORT POL 


This subroutine, from the text by Conte and de Boor, 
{2] forms a fifth order least squares polynomial through a 


set of data points, representing values of Q(k ) and k=, 

iG is 
each with a weight W. fhe resulting polynomial, evaluated 
by the subroutine ORTGET, finds the first estimate of k 


re 
for a given value of Q. ORTGET returns a value of the 


polynomial and its derivative in terms of both phase _ speed 
(CPV) and horizontal wave number (KM), as well as the 
derivative of the k with respect to Q (DKRDQ). 

cr 


5- EIGFUN 





Given a value of k and a depth, Zz, (entered by 
Een 


depth grid index (IDZR) and the depth below the grid depth 


(DELZR)), the subroutine EIGFUN calculates the unnormalized 
eigenfunction U(zZ).~ Upon the first call to EIGFUN after the 


determination of the eigenvalue, the WKB coefficients, a or 
O 


S (CO1) and b or d {CO2) are calculated for each depth 
O O O 
grid as a function of the upper and lower turning point 


phase angles (THU and THL). The WKB and Airy validity 
factors are evaluated; if either is less than the parameters 
KALIM or KWLIM, the Airy or WKB solution is calculated 
immediately. If neither the WKB or Airy conditions is 
satisfied, a search is made both above and below z, to find 
the first depths with valid WKB or Airy solutions. At these 
two depths U and JU are calculated; then the values are 


passed to the subroutine SPLICE. 
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6. SPLICE 


Starting with values of U and U* (U and UP), SPLICE 


extrapolates between 2Z and a (DZP) using a finite 
a 


difference version of Eq. (2.13), 


Zoe 
Gee eeeriz) aie 0 0 
T+ 1 Z 1 y= 
First an extrapolation starts from the lower depth and 
iterates in 101 or less steps towards the upper depth. This 
finite difference scheme is guite unstable in the shadow 


zone region and the extrapolation function generated, U1), 
1 


may exponentially grow with depth. In order to force U_ to 
i 


be continuous at Ze 022 ais now multiplied by a ramp 
a 


function forning, 


meee YCtset (0 = USLIX) (2 -z ) / (2 -Z +). 
a 1 D 1 1a boa 


This finite difference scheme is repeated in the opposite 


direction (fron Zs to a to form an extrapolation function 
_ This in turn is multiplied by a second ramp function 
yielding the extrapolation function ae which converges to 
mz} There are now two continuous extrapolation functions 
a" and re with ace representing a more accurate 
extrapolation near Zi and a representing a more accurate 
extrapolation near aise An average of the two extrapolators 


weighted by the relative proximity of the desired depth to 


Z and ze is used for the final estimate of the 
a 
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eigenfunction, 
Peo | (2aaey USS2e se (7 —-2z) US4)) 7 (z -z). 
n b 1 a 1 a b 


Due to the ramp function used, this extrapolation scheme may 


give a discontinuous derivative at the end points, 2z and 
a 


Za Since the derivative of U(z) with depth plays no part 


in the calculation of the transmission loss, and the errors 
incurred in U(z) are slight, this effect can be expected to 


have little if any effect on the subsequent calculations. 
fe 2LOSs 


The subroutine PLOSS computes the transmission loss 
With range from RSTART to RSTOP in increments of DELR, and 
displays the resultsS on a printer graph or a CALCOMP plot. 
The transmission loss is caiculated in terms of both the 
coherent and incoherent transmission loss, and an optional 
averaged transmission loss. The averaged transmission loss 


is calculated from a running average of the square of the 


2 
coherent pressure P multiplied by a range weighting 


function (or window). Two windows are provided, a Gaussian 


curve (IFILT=1) and a straight running average (IFILT=2). 
When PLOSS is called, various initial values are set 
and the filter function is normalized. The averaging is 


2 
performed with a stack which is filled with P values. At 
each range step the values in the stack drop one place; the 


- . 
value of P is calculated for the range r plus the averaging 


width (RAVG) and is placed at the end of the stack. The 
coherent transmission loss is calculated from the middle 


value of the stack. 
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For multiple profile plots the transmission loss 
calculations stop at the range equal to the end of the 
profile sector (REND), and a variable is set indicating the 
Subroutine was left "on" ($0N). After processing of the 
next profile, the transmission loss calculations continue 


without interruption to the output. 


8. MATCH 


This subroutine fills two vectors (ICROSS,IREX) 
which provide a cross-reference table for the matching of 
modes between profiles. fhe matching is carried out in two 
steps: 

s Each mode in the new profile is selected in order of 
increasing phase speed (MORDR2). For each mode, the modes 
in the previous profile are searched (also in order of 
increasing phase speed) for a mode with the same relative 
position in the profile ({IBRN1,IBRN2) and with the same 
mode number (MODEL,MODEN). 

a iIf, after the first stev, any modes remain unmatched 
a transition between channels has taken place. In this 
case the second step begins by finding the extent of the 
gap of unmatched modes. The gap begins with the first 
unmatched mode in the new profile, and ends with the next 
matched mode in either the same channei or the nain 
channel. The corresponding gap is aiso found for the 
previous profile, and the number of missing modes in each 
channel is tabulated. The modes are then cross-referenced 
(ICROSS) in terms of increasing phase speed. As long as 
unmatched modes remain in more than one sub-channel, the 
source mode energy is summed (USUM), later to be divided 
evenly between the same nodes. 

This step is repeated for all such gaps. Finally, any modes 
which still remain unmatched are cross-referenced to mode 


200, which is used aS a dummy, or null mode. 


128 





9. INTEGR 


The subroutine INTEGR integrates either the vertical 
real wave number or vertical imaginary wave number with 
depth. With the assumed nature of the profile this 
integration is accomplished in each layer with the analytic 


formula, 
ance 3 3 
kK dz = 2 {i(k (z.) - kK (z, 1) WH /y- 
Z Z 3 Zz i Z 


10. QPIND 


The subroutine QFIND computes the value of Q at a 


particular k for a family between the depth grids MINDEP 
Fe 


and MAXDEP. This 1s accomplished by integrating the 


vertical wave number across the ensonified depths (INTEGR), 
then calling the routines for the surface and  botton 


boundary conditions (SBC and BBC). 
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VIII. RESULTS AND CONCLUSION 


Results form the QMODE model were compared with 
analytic solutions, the results from some other models, and 


with some transmission loss experimental data. 


A. ANALYTIC TESTS 


The model was tested with two profiles which yield 


analytic values for the eigenvalues and eigenfunctions. 


1. Isovelocity Channel 


The first profile is that of an isovelocity sound 
channel With rigid bottom, for which the eigenvalues are 
2 Z Z 

K =k - ((2n-1) /(25)] , (8.1) 

r,n 
where H is the depth of the channel. The profile hada 
sound speed of 1500 n/fsec and a depth of 4000 meters. The 
first 150 eigenvalues were computed for a frequency of 100 
Hz by both the QMODE and NORMO1 (a normal anode finite 
difference program by the author {1]) programs. The results 
are shown in Fig 23, where the error between the analytic 
and computed eigenvalues, scaled to the intermode spacing, 
is plotted as function of mode number. A curve is also 
included which corresponds to six significant figure 
agreement between the computed eigenvalues and the analytic 
eigenvalues. For higher nodes, as the vertical grid spacing 
becomes an appreciable fraction of the vertical wavelength, 


the finite difference technique becomes more inaccurate. 
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2. Constant Gradient Channel 

The second profile is one with a single c gradient 
(Fig 29). For this profile the analytic eigenfunctions are 
a linear combination of the Airy functions, Ai and Bi. [In 
the shadow zone Bi increases exponentially with depth, for 
modes with lower turning point well above the botton, and 
the coefficient of the Bi term must be zero. Thus the 
eigenvalues are those values of the horizontal wave number 


for which the zeros of the Airy function occur at the 


surface. This occurs when (using the notation of chapter 
II), 
Ai(-%) = 0, (8.2) 
oe, 
2 23 2 
K = ¥ Va ce, (8. 3) 
nea n O 


where y are the negative zeros of the Airy function and k 
n O 


is the wave number at the surface. The differences between 
the analytic eigenvalues and the computed eigenvalues are 
shown in Fig 30 as a function of mode number. For the lower 
modes the finite difference technique and the WKB model have 
comparable accuracy. Beginning with mode 12 the QMODE (WKB) 
eigenvalues become accurate to six significant figures, 


asymptoticaily approaching the analytic solution. 


Be. AESD WORKSHOP TESTS 


The QMODE results for three profiles used in the AESD 
Workshop on Non-Ray-Tracing Techniques were examined [2]. 
The three profiles correspond to a Pacific deep water 
profile (AESD 1B), a N.E. Atlantic profile with multiple 


sound channeis (AESD 2), and a Shallow water profile 
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(AESD 3). The sound speed profiles for each are shown in 
pig 31, 32, and 33. 


1. Comparison of Modes 


= SSS SSS SSS SS 


From two of these profiles, a number of 
eigenfunctions were calculated by both the QMODE and NORMO1 
programs and plotted. All five eigenfunctions for the 
Shallow water case are shown in Fig 34. fhe normalized 
eigenfunctions generated by both programs are superimposed, 
with the NORMO1 modes in solid lines, and the QMODE modes in 
dashed lines. Modes 11 through 20 of the multiple channel 
profile, AESD 2, are shown in Fig 35 and 36. fhese modes 
occurred at phase speeds near the edge of the barrier 


separating tne two channels. 
2. Transmission Loss 


Transmission loss was calculated by QMODE for the 
AESD profiies. The QMODE results are compared with those 
for the FACT (modified ray theory), the NOL (finite 
difference normal mode), and the parabolic equation models. 
A description of the NOL modei is found in Ref. {3]. The 
central processing unit (CPU) times for each QMODE run are 
given in Tables 1 through 3. The transmission loss curves 
are given beginning with Pics 37. The QMODE transmission 
loss plots may show up to three curves: the fully coherent 
transmission loss, the incoherent transmission loss, and the 
averaged transmission loss. The QMODE fresultS agree 
generally, although not necessarily in detail, with those 
for the finite difference and parabolic equation models. 
The FACT results are considerably smoothed and do not 
reflect the sharp interference effectS apparent in the 
normal mode and parabolic equation nodels. The FACT model 
incorporates acoustic energy with higher grazing angles at 
the bottom (greater than critical angle) than normal mode 
techniques. Thus the QMODE results show greater 


transmission loss in the shadow zones, especially for 
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profile AESD 1B. 

If RAVG, the range over which the QMODE transmission 
loss is averaged, is increased the results are smoothed and 
begin to resemble the FACT results. This is illustrated in 
Fig 49 through Fig 52, where the averaging range for the 
AESD 1B transmission loss is increased from 2.5 nm to 10 


nM. 


C. COMPARISONS WITH DATA 


The transmission loss predicted by QMODE was compared 


with the results from two sets of acoustic experiments. 


1. PARKA IIA Data 


The first set of comparisons was with a series of 
experiments in the North Pacific (PARKA IIA) [4]. The 
experiment measured transmission Loss over different paths 
for various combinations of frequency, source, and frecelver 
depth. The research platform FLIP, equipped with suspended 
hydrophones, acted as the receiver site. Explosive charges 
were dropped along radial bearings from the receiver and the 
resulting transmission loss was recorded. Four runs were 
compared, two to ranges of 500 nm, one to 1500 nm, and one 
to 2000 an. 


ae Run 1 


In this run the transmission loss at 100 Hz 
between a deeo source (500 ft) and a sound channel axis 
receiver (2563 ft) was measured out to a range of 500 nna. 
The results are shown in Fig 53. The QMODE results agree 
well in the convergence zones and at longer ranges, but snow 
greater transmission loss in the shadow zones between the 
first few convergence zones. This can be attributed to the 
fact that QMODES was run without including any bottom bounce 


energy (Since detailed bottom information was not readily 
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available). 
bee Aun 2 


This example measured transmission Loss at 100 
HZ between a shallow source (60 £t) and a shallow (300 ft) 
receiver. The data were collected in two runs, one outbound 
from the receiver for a distance of 500 nm, and then a run 
back towards the receiver. The results for the outbound and 
inbound runs are shown in Pig 54 and 55. Again the QMODE 
results agree with the data at the longer ranges and in the 
convergence zones. The effect of omitting bottom bounce 


energy is still apparent. 
CGeaekUn 3 


The third PARKA example measured the 
transmission loss at 100 Hz between a deep source (600 ft) 
and a sound channel axis receiver out to a distance of 1500 
nm. The environmental input for the QMODE computer run was 
based upon the Pleet Numerical Weather Central (FNWC) 
climatology for the appropriate month. This climatology is 
formed by averaging the temperature and salinity fields for 
historical data. The computer run used a profile based on 
Climatology for every 100 nm along the track. The results 
are shown in Fig 56. For this northerly track the sound 
channel siowly rises with range and the snallow source 


becomes increasingly coupled to the axis receiver. 
de Run 4 


The last example from the PARKA experiments was 
along an east-west track for a deep (800 ft) source and axis 
receiver out to a distance of 2000 nm. Along this track the 
sound channel axis depth remained relatively constant. The 
environmental data for the computer run was based upon tne 
FNWC climatology with the temperature values corrected at 
Shallow depths for observed AXBI (Airborne Expendable 
Bathythermograph) data collected at the time of the 
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experiment [5]. The results are shown in Fig 57. For the 
first 1000 miles the QMODE results agree well with the 
Observed results. After this range the QMODE model fails to 
predict the observed large variations in the transmission 
loss. The FNWC propagation loss calculations for the same 
track using the RP-70 model (ray-~-trace) also failed to show 
the observed variations [4]. It is assumed that the wide 
variation of transmission loss is due to differences between 


the assumed and actual bathymetry along the track. ° 
2. Antigua to Grand Banks 


The second acoustic experiment for which comparisons 
were made measured the transmission loss along a great 
circle track from a location near Antigua, West Indies to 
the Grand Banks. The results of this experiment were 
described in an article by Guthrie, Fitzgerald, Nutile, and 
Shaffer [6]. The author was present at the receiver site in 
connection with military duties. Two acoustic CW projectors 
at frequencies of 13.89 and 111.1 Hz were towed at depths 
of 104 and 21 meters. The received signal was recorded 
through a 1100-meter deep, bottom-moored hydrophone off 
Antigua. The QMODE results show the effects of variable 
oceanographic conditions along the track. For each source 
frequency two sets of calculations were made, one With a 
bilinear profile (Fig 58), and a set of profiles based upon 
FNWC climatology at 100-aom intervals. For the low frequency 
source, Fig 59, 60, and 61 show the results for the bilinear 
profile and the climatology profiles, as well as the 
observed transmission loss. The QMODE results show 
well-defined convergence zones along the entire length of 
the track, while the observed data show the zones to lose 
definition after about 1800 ka. 

Pigures 62, 63 and 64 show the corresponding results 
for the high-frequency source. Again, QMODE predicted 
well-defined convergence zones to greater ranges than was 


observed. The appearance of well defined modes at greater 
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ranges suggests that more scattering of acoustic energy 
among the modes took place than was accounted for in the 
computer calculations. This scattering can be due to either 
inhomgeneities in the ocean or the non-adiabatic coupling 
among modes of different profiles. Guthrie, et al, observed 
that as the track extends to the north the sound channel 
becomes shallower and more acoustic energy is coupled fron 
the shallow sources into the deep sound channel. This 
change in the nature of the sound speed profile is observed 
in the QMODE results, where the incoherent transmission loss 


decreases as the northern latitudes are reached. 


D. EFFECT OF THE BOTTOM ON DEEP SOUND CHANNEL PROPAGATION 


In the PARKA run 4, lack of information on the botton 
appeared to have an adverse effect on the computer results. 
Due to a lack of data, the computer runs were made 
considering only refracted energy (this approxigation is 
often used when dealing with such long range transmission). 
Along tracks where the modes have little interaction with 
the bottom, this limitation has little effect. This usually 
occurs in the Northern Hemisphere for northerly tracks 
where, typically, the sound channel becomes shallower with 
increasing latitude. On an east-west track (as in PARKA run 
4), and on soutnerly tracks (where the sound channel becomes 
deeper with lower latitudes), the modes tend to interact 
more with the bottom. For such tracks the calculation of 
deep sound channel transmission loss becomes difficult 
without detailed knowledge of the bottom characteristics. 
Thus to predict low-frequency sound propagation reliably 
over long distances in the ocean, detailed botton 
information is required in a form suitable for the 


mathematical model being used. 
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Ee. CONCLUSIONS 


QMODE possesses the potential of a relatively fast 
method for the calculation of long range low-frequency sound 
transmission in the oceans. The computational speed results 
from the normal mode method of computing transmission loss 
and the method in which the eigenvalues and eigenfunctions 
(modes) are determined. | 

The theoretical basis for the eigenvalues and the 
eigenfunctions is the assumed WKB solution. The first 
estimate for the eigenvalue is obtained by using a 
polynomial estimator through data generated by the WKB 
characteristic equation. When a precise eigenvalue is 
required, this estimator allows the use of significantly 
fewer iterations. For the calculation of averaged or 
incoherent transmission loss, the estimator alone provides 
an accurate enough approximation of the eigenvalue. Since 
the eigenfunction is an assumed WKB function, a very precise 
elgevalue iS not freguired in order to generate a good 
eigenfunction. This contrasts with conditions in the finite 
difference technique, where a precise eigenvalue is required 
due to instabilities in the finite difference equation in 
the shadow zones. 

The second characteristic of the model which aliows 
rapid computational speed 1S a property common to all normal 
mode models: the accuracy of the transmission loss 
calculations does not depend upon the size of the range step 
between calculations. fo calculate the transmission loss at 
some range, the following is required: 

a The modes be properly defined at the receiver; 

s The coupling between modes be evaluated in a range 
dependent environment, 

es The mode attenuation be evaluated; 


=» The modes be defined at the receiver. 


Thus, the steps required in this model are dependent 
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upon the number of intervening profiles and not the range 
between source and receiver. In contrast, cay-tracing 
methods must calculate ray propagation over small grid steps 
along the entire path. A limitation of the mode method is 
that the number of modes, hence computation time, is 
proportional to the frequency considered. 

There are a number of applications for which the model 
is ill-suited. The WKB method should not be expected to 
yield accurate results in situations in which only a few 
modeS are present, such as in a sound channel near the 
cut~off frequency. As has been seen, the few modes near the 
edge of a barrier which separates multiple sound channels 
Should not be expected to be exact. In the computation of 
transmission loss, this should not be a limitation if more 
than a few modes are present. A weakness of the normal mode 
technique is that rapid variations of sound speed and botton 
depth with range are not easily accommodated. A high 
relief, sharp bottom becomes extremely difficult to model; 
this is true as well for rapid changes in the sound speed 
profile. 

A number of improvements concerning both the theory and 
the actual program could be made to the model as presently 
configured. Some method of eStimating the effects of a 
steep~sloped bottom upon the distribution of energy in the 
modes would improve the model's ability to predict 
transmission loss over variable bathymetry. Such an 
estimate would physically correspond to the non-adiabatic 
coupling between modes; or in terms of rays, the ray 
deflection after striking a sloping botton. Secondly, an 
estimate of the scattering of energy in the modes due to 
both inhomogeneities in the environment and to the 
non~adiabatic transition between profiles would be 
desirable. A number of programming improvements could be 
made, including more efficient coding, space for storing a 
larger number of modes (presently 200), and multiple 


source/receiver deptas. 
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FIGURE 28 - Eigenvalue Test: Isovelocity Channel 
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